为啥这个 SIMD 乘法不比非 SIMD 乘法快?
Posted
技术标签:
【中文标题】为啥这个 SIMD 乘法不比非 SIMD 乘法快?【英文标题】:Why is this SIMD multiplication not faster than non-SIMD multiplication?为什么这个 SIMD 乘法不比非 SIMD 乘法快? 【发布时间】:2017-08-15 08:16:30 【问题描述】:假设我们有一个函数将两个数组相乘,每个数组有 1000000 个双精度数。在 C/C++ 中,函数如下所示:
void mul_c(double* a, double* b)
for (int i = 0; i != 1000000; ++i)
a[i] = a[i] * b[i];
编译器使用-O2
生成以下程序集:
mul_c(double*, double*):
xor eax, eax
.L2:
movsd xmm0, QWORD PTR [rdi+rax]
mulsd xmm0, QWORD PTR [rsi+rax]
movsd QWORD PTR [rdi+rax], xmm0
add rax, 8
cmp rax, 8000000
jne .L2
rep ret
从上面的程序集中,编译器似乎使用了 SIMD 指令,但它每次迭代只乘以一倍。所以我决定改为在内联汇编中编写相同的函数,在那里我充分利用xmm0
寄存器并一次将两个双精度数相乘:
void mul_asm(double* a, double* b)
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"0: \n\t"
"movupd xmm0, xmmword ptr [rdi+rax] \n\t"
"mulpd xmm0, xmmword ptr [rsi+rax] \n\t"
"movupd xmmword ptr [rdi+rax], xmm0 \n\t"
"add rax, 16 \n\t"
"cmp rax, 8000000 \n\t"
"jne 0b \n\t"
".att_syntax noprefix \n\t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
分别测量这两个函数的执行时间后,似乎它们都需要 1 ms 才能完成:
> gcc -O2 main.cpp
> ./a.out < input
mul_c: 1 ms
mul_asm: 1 ms
[a lot of doubles...]
我希望 SIMD 实现的速度至少是乘法/内存指令数量的一半(0 ms)的两倍。
所以我的问题是:当 SIMD 实现只执行一半的乘法/内存指令时,为什么 SIMD 实现不比普通 C/C++ 实现快?
这是完整的程序:
#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
void mul_c(double* a, double* b)
for (int i = 0; i != 1000000; ++i)
a[i] = a[i] * b[i];
void mul_asm(double* a, double* b)
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"0: \n\t"
"movupd xmm0, xmmword ptr [rdi+rax] \n\t"
"mulpd xmm0, xmmword ptr [rsi+rax] \n\t"
"movupd xmmword ptr [rdi+rax], xmm0 \n\t"
"add rax, 16 \n\t"
"cmp rax, 8000000 \n\t"
"jne 0b \n\t"
".att_syntax noprefix \n\t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
int main()
struct timeval t1;
struct timeval t2;
unsigned long long time;
double* a = (double*)malloc(sizeof(double) * 1000000);
double* b = (double*)malloc(sizeof(double) * 1000000);
double* c = (double*)malloc(sizeof(double) * 1000000);
for (int i = 0; i != 1000000; ++i)
double v;
scanf("%lf", &v);
a[i] = v;
b[i] = v;
c[i] = v;
gettimeofday(&t1, NULL);
mul_c(a, b);
gettimeofday(&t2, NULL);
time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
printf("mul_c: %llu ms\n", time);
gettimeofday(&t1, NULL);
mul_asm(b, c);
gettimeofday(&t2, NULL);
time = 1000 * (t2.tv_sec - t1.tv_sec) + (t2.tv_usec - t1.tv_usec) / 1000;
printf("mul_asm: %llu ms\n\n", time);
for (int i = 0; i != 1000000; ++i)
printf("%lf\t\t\t%lf\n", a[i], b[i]);
return 0;
我还尝试利用所有xmm
寄存器(0-7)并删除指令依赖以获得更好的并行计算:
void mul_asm(double* a, double* b)
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"0: \n\t"
"movupd xmm0, xmmword ptr [rdi+rax] \n\t"
"movupd xmm1, xmmword ptr [rdi+rax+16] \n\t"
"movupd xmm2, xmmword ptr [rdi+rax+32] \n\t"
"movupd xmm3, xmmword ptr [rdi+rax+48] \n\t"
"movupd xmm4, xmmword ptr [rdi+rax+64] \n\t"
"movupd xmm5, xmmword ptr [rdi+rax+80] \n\t"
"movupd xmm6, xmmword ptr [rdi+rax+96] \n\t"
"movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
"mulpd xmm0, xmmword ptr [rsi+rax] \n\t"
"mulpd xmm1, xmmword ptr [rsi+rax+16] \n\t"
"mulpd xmm2, xmmword ptr [rsi+rax+32] \n\t"
"mulpd xmm3, xmmword ptr [rsi+rax+48] \n\t"
"mulpd xmm4, xmmword ptr [rsi+rax+64] \n\t"
"mulpd xmm5, xmmword ptr [rsi+rax+80] \n\t"
"mulpd xmm6, xmmword ptr [rsi+rax+96] \n\t"
"mulpd xmm7, xmmword ptr [rsi+rax+112] \n\t"
"movupd xmmword ptr [rdi+rax], xmm0 \n\t"
"movupd xmmword ptr [rdi+rax+16], xmm1 \n\t"
"movupd xmmword ptr [rdi+rax+32], xmm2 \n\t"
"movupd xmmword ptr [rdi+rax+48], xmm3 \n\t"
"movupd xmmword ptr [rdi+rax+64], xmm4 \n\t"
"movupd xmmword ptr [rdi+rax+80], xmm5 \n\t"
"movupd xmmword ptr [rdi+rax+96], xmm6 \n\t"
"movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
"add rax, 128 \n\t"
"cmp rax, 8000000 \n\t"
"jne 0b \n\t"
".att_syntax noprefix \n\t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
但它仍然以 1 毫秒的速度运行,与普通 C/C++ 实现的速度相同。
更新
正如 answers/cmets 所建议的,我已经实现了另一种测量执行时间的方法:
#include <stdio.h>
#include <stdlib.h>
void mul_c(double* a, double* b)
for (int i = 0; i != 1000000; ++i)
a[i] = a[i] * b[i];
void mul_asm(double* a, double* b)
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"0: \n\t"
"movupd xmm0, xmmword ptr [rdi+rax] \n\t"
"mulpd xmm0, xmmword ptr [rsi+rax] \n\t"
"movupd xmmword ptr [rdi+rax], xmm0 \n\t"
"add rax, 16 \n\t"
"cmp rax, 8000000 \n\t"
"jne 0b \n\t"
".att_syntax noprefix \n\t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
void mul_asm2(double* a, double* b)
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"0: \n\t"
"movupd xmm0, xmmword ptr [rdi+rax] \n\t"
"movupd xmm1, xmmword ptr [rdi+rax+16] \n\t"
"movupd xmm2, xmmword ptr [rdi+rax+32] \n\t"
"movupd xmm3, xmmword ptr [rdi+rax+48] \n\t"
"movupd xmm4, xmmword ptr [rdi+rax+64] \n\t"
"movupd xmm5, xmmword ptr [rdi+rax+80] \n\t"
"movupd xmm6, xmmword ptr [rdi+rax+96] \n\t"
"movupd xmm7, xmmword ptr [rdi+rax+112] \n\t"
"mulpd xmm0, xmmword ptr [rsi+rax] \n\t"
"mulpd xmm1, xmmword ptr [rsi+rax+16] \n\t"
"mulpd xmm2, xmmword ptr [rsi+rax+32] \n\t"
"mulpd xmm3, xmmword ptr [rsi+rax+48] \n\t"
"mulpd xmm4, xmmword ptr [rsi+rax+64] \n\t"
"mulpd xmm5, xmmword ptr [rsi+rax+80] \n\t"
"mulpd xmm6, xmmword ptr [rsi+rax+96] \n\t"
"mulpd xmm7, xmmword ptr [rsi+rax+112] \n\t"
"movupd xmmword ptr [rdi+rax], xmm0 \n\t"
"movupd xmmword ptr [rdi+rax+16], xmm1 \n\t"
"movupd xmmword ptr [rdi+rax+32], xmm2 \n\t"
"movupd xmmword ptr [rdi+rax+48], xmm3 \n\t"
"movupd xmmword ptr [rdi+rax+64], xmm4 \n\t"
"movupd xmmword ptr [rdi+rax+80], xmm5 \n\t"
"movupd xmmword ptr [rdi+rax+96], xmm6 \n\t"
"movupd xmmword ptr [rdi+rax+112], xmm7 \n\t"
"add rax, 128 \n\t"
"cmp rax, 8000000 \n\t"
"jne 0b \n\t"
".att_syntax noprefix \n\t"
:
: "D" (a), "S" (b)
: "memory", "cc"
);
unsigned long timestamp()
unsigned long a;
asm volatile
(
".intel_syntax noprefix \n\t"
"xor rax, rax \n\t"
"xor rdx, rdx \n\t"
"RDTSCP \n\t"
"shl rdx, 32 \n\t"
"or rax, rdx \n\t"
".att_syntax noprefix \n\t"
: "=a" (a)
:
: "memory", "cc"
);
return a;
int main()
unsigned long t1;
unsigned long t2;
double* a;
double* b;
a = (double*)malloc(sizeof(double) * 1000000);
b = (double*)malloc(sizeof(double) * 1000000);
for (int i = 0; i != 1000000; ++i)
double v;
scanf("%lf", &v);
a[i] = v;
b[i] = v;
t1 = timestamp();
mul_c(a, b);
//mul_asm(a, b);
//mul_asm2(a, b);
t2 = timestamp();
printf("mul_c: %lu cycles\n\n", t2 - t1);
for (int i = 0; i != 1000000; ++i)
printf("%lf\t\t\t%lf\n", a[i], b[i]);
return 0;
当我用这个测量值运行程序时,我得到了这个结果:
mul_c: ~2163971628 cycles
mul_asm: ~2532045184 cycles
mul_asm2: ~5230488 cycles <-- what???
这里有两点值得注意,首先,周期数变化很大,我认为这是因为操作系统允许其他进程在其间运行。有什么方法可以防止这种情况发生,或者只在我的程序执行时计算周期?此外,mul_asm2
产生的输出与其他两个相同,但速度要快得多,如何?
我在我的系统上尝试了 Z boson 的程序以及我的 2 个实现,得到了以下结果:
> g++ -O2 -fopenmp main.cpp
> ./a.out
mul time 1.33, 18.08 GB/s
mul_SSE time 1.13, 21.24 GB/s
mul_SSE_NT time 1.51, 15.88 GB/s
mul_SSE_OMP time 0.79, 30.28 GB/s
mul_SSE_v2 time 1.12, 21.49 GB/s
mul_v2 time 1.26, 18.99 GB/s
mul_asm time 1.12, 21.50 GB/s
mul_asm2 time 1.09, 22.08 GB/s
【问题讨论】:
您的时间计算对于这种基准测试来说不够精确。尝试使用Google Benchmark library 运行代码,看看你会发现什么。 您需要更多的循环迭代才能更好地测量它,使用高分辨率计时器或使用 RDTSC/RDTSCP。你有 1 毫秒是噪音。 比如你可能被内存瓶颈了。 另外使用 -O3,你将拥有 C 版本的mulpd xmm0, XMMWORD PTR [rcx+rax]
。
你这里的内存绝对是瓶颈。
【参考方案1】:
以前的基准测试有 a major bug in the timing function I used。这严重低估了没有矢量化和其他测量的带宽。此外,还有另一个问题是高估了阵列上已读取但未写入的带宽due to COW。最后,我使用的最大带宽不正确。我已通过更正更新了我的答案,并在此答案的末尾留下了旧答案。
您的操作受内存带宽限制。这意味着 CPU 大部分时间都在等待缓慢的内存读取和写入。可以在这里找到一个很好的解释:Why vectorizing the loop does not have performance improvement。
但是,我不得不稍微不同意那个答案中的一个陈述。
因此,无论如何优化(矢量化、展开等),它都不会变得更快。
事实上,即使在内存带宽受限的操作中,矢量化、展开和多线程也可以显着增加带宽。原因是难以获得最大的内存带宽。可以在这里找到一个很好的解释:https://***.com/a/25187492/2542702。
我的其余答案将展示矢量化和多线程如何接近最大内存带宽。
我的测试系统:Ubuntu 16.10,Skylake (i7-6700HQ@2.60GHz),32GB RAM,双通道 DDR4@2400 GHz。我系统的最大带宽为 38.4 GB/s。
根据下面的代码,我生成了下表。我使用 OMP_NUM_THREADS 例如设置线程数export OMP_NUM_THREADS=4
。效率为bandwidth/max_bandwidth
。
-O2 -march=native -fopenmp
Threads Efficiency
1 59.2%
2 76.6%
4 74.3%
8 70.7%
-O2 -march=native -fopenmp -funroll-loops
1 55.8%
2 76.5%
4 72.1%
8 72.2%
-O3 -march=native -fopenmp
1 63.9%
2 74.6%
4 63.9%
8 63.2%
-O3 -march=native -fopenmp -mprefer-avx128
1 67.8%
2 76.0%
4 63.9%
8 63.2%
-O3 -march=native -fopenmp -mprefer-avx128 -funroll-loops
1 68.8%
2 73.9%
4 69.0%
8 66.8%
由于测量中存在不确定性,经过多次迭代运行,我得出以下结论:
单线程标量操作获得超过 50% 的带宽。 两个线程标量操作获得最高带宽。 单线程向量操作比单线程标量操作更快。 单线程 SSE 操作比单线程 AVX 操作更快。 展开没有帮助。 展开单线程操作比不展开要慢。 线程数多于内核数(超线程)会降低带宽。提供最佳带宽的解决方案是具有两个线程的标量操作。
我用来基准测试的代码:
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <omp.h>
#define N 10000000
#define R 100
void mul(double *a, double *b)
#pragma omp parallel for
for (int i = 0; i<N; i++) a[i] *= b[i];
int main()
double maxbw = 2.4*2*8; // 2.4GHz * 2-channels * 64-bits * 1-byte/8-bits
double mem = 3*sizeof(double)*N*R*1E-9; // GB
double *a = (double*)malloc(sizeof *a * N);
double *b = (double*)malloc(sizeof *b * N);
//due to copy-on-write b must be initialized to get the correct bandwidth
//also, GCC will convert malloc + memset(0) to calloc so use memset(1)
memset(b, 1, sizeof *b * N);
double dtime = -omp_get_wtime();
for(int i=0; i<R; i++) mul(a,b);
dtime += omp_get_wtime();
printf("%.2f s, %.1f GB/s, %.1f%%\n", dtime, mem/dtime, 100*mem/dtime/maxbw);
free(a), free(b);
有计时错误的旧解决方案
内联汇编的现代解决方案是使用内部函数。仍然存在需要内联汇编的情况,但这不是其中之一。
内联汇编方法的一个内在解决方案很简单:
void mul_SSE(double* a, double* b)
for (int i = 0; i<N/2; i++)
_mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
让我定义一些测试代码
#include <x86intrin.h>
#include <string.h>
#include <stdio.h>
#include <x86intrin.h>
#include <omp.h>
#define N 1000000
#define R 1000
typedef __attribute__(( aligned(32))) double aligned_double;
void (*fp)(aligned_double *a, aligned_double *b);
void mul(aligned_double* __restrict a, aligned_double* __restrict b)
for (int i = 0; i<N; i++) a[i] *= b[i];
void mul_SSE(double* a, double* b)
for (int i = 0; i<N/2; i++) _mm_store_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
void mul_SSE_NT(double* a, double* b)
for (int i = 0; i<N/2; i++) _mm_stream_pd(&a[2*i], _mm_mul_pd(_mm_load_pd(&a[2*i]),_mm_load_pd(&b[2*i])));
void mul_SSE_OMP(double* a, double* b)
#pragma omp parallel for
for (int i = 0; i<N; i++) a[i] *= b[i];
void test(aligned_double *a, aligned_double *b, const char *name)
double dtime;
const double mem = 3*sizeof(double)*N*R/1024/1024/1024;
const double maxbw = 34.1;
dtime = -omp_get_wtime();
for(int i=0; i<R; i++) fp(a,b);
dtime += omp_get_wtime();
printf("%s \t time %.2f s, %.1f GB/s, efficency %.1f%%\n", name, dtime, mem/dtime, 100*mem/dtime/maxbw);
int main()
double *a = (double*)_mm_malloc(sizeof *a * N, 32);
double *b = (double*)_mm_malloc(sizeof *b * N, 32);
//b must be initialized to get the correct bandwidth!!!
memset(a, 1, sizeof *a * N);
memset(b, 1, sizeof *a * N);
fp = mul, test(a,b, "mul ");
fp = mul_SSE, test(a,b, "mul_SSE ");
fp = mul_SSE_NT, test(a,b, "mul_SSE_NT ");
fp = mul_SSE_OMP, test(a,b, "mul_SSE_OMP");
_mm_free(a), _mm_free(b);
现在是第一个测试
g++ -O2 -fopenmp test.cpp
./a.out
mul time 1.67 s, 13.1 GB/s, efficiency 38.5%
mul_SSE time 1.00 s, 21.9 GB/s, efficiency 64.3%
mul_SSE_NT time 1.05 s, 20.9 GB/s, efficiency 61.4%
mul_SSE_OMP time 0.74 s, 29.7 GB/s, efficiency 87.0%
所以-O2
不向量化循环,我们看到内在 SSE 版本比普通 C 解决方案 mul
快得多。 efficiency = bandwith_measured/max_bandwidth
我的系统的最大值为 34.1 GB/s。
第二次测试
g++ -O3 -fopenmp test.cpp
./a.out
mul time 1.05 s, 20.9 GB/s, efficiency 61.2%
mul_SSE time 0.99 s, 22.3 GB/s, efficiency 65.3%
mul_SSE_NT time 1.01 s, 21.7 GB/s, efficiency 63.7%
mul_SSE_OMP time 0.68 s, 32.5 GB/s, efficiency 95.2%
-O3
对循环进行矢量化处理,而内部函数基本上没有任何优势。
第三次测试
g++ -O3 -fopenmp -funroll-loops test.cpp
./a.out
mul time 0.85 s, 25.9 GB/s, efficency 76.1%
mul_SSE time 0.84 s, 26.2 GB/s, efficency 76.7%
mul_SSE_NT time 1.06 s, 20.8 GB/s, efficency 61.0%
mul_SSE_OMP time 0.76 s, 29.0 GB/s, efficency 85.0%
使用-funroll-loops
GCC 将循环展开八次,除了非临时存储解决方案和 OpenMP 解决方案没有真正优势外,我们看到了显着的改进。
在展开循环之前,mul
和 -O3
的程序集是
xor eax, eax
.L2:
movupd xmm0, XMMWORD PTR [rsi+rax]
mulpd xmm0, XMMWORD PTR [rdi+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
add rax, 16
cmp rax, 8000000
jne .L2
rep ret
对于-O3 -funroll-loops
,mul
的程序集是:
xor eax, eax
.L2:
movupd xmm0, XMMWORD PTR [rsi+rax]
movupd xmm1, XMMWORD PTR [rsi+16+rax]
mulpd xmm0, XMMWORD PTR [rdi+rax]
movupd xmm2, XMMWORD PTR [rsi+32+rax]
mulpd xmm1, XMMWORD PTR [rdi+16+rax]
movupd xmm3, XMMWORD PTR [rsi+48+rax]
mulpd xmm2, XMMWORD PTR [rdi+32+rax]
movupd xmm4, XMMWORD PTR [rsi+64+rax]
mulpd xmm3, XMMWORD PTR [rdi+48+rax]
movupd xmm5, XMMWORD PTR [rsi+80+rax]
mulpd xmm4, XMMWORD PTR [rdi+64+rax]
movupd xmm6, XMMWORD PTR [rsi+96+rax]
mulpd xmm5, XMMWORD PTR [rdi+80+rax]
movupd xmm7, XMMWORD PTR [rsi+112+rax]
mulpd xmm6, XMMWORD PTR [rdi+96+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
mulpd xmm7, XMMWORD PTR [rdi+112+rax]
movaps XMMWORD PTR [rdi+16+rax], xmm1
movaps XMMWORD PTR [rdi+32+rax], xmm2
movaps XMMWORD PTR [rdi+48+rax], xmm3
movaps XMMWORD PTR [rdi+64+rax], xmm4
movaps XMMWORD PTR [rdi+80+rax], xmm5
movaps XMMWORD PTR [rdi+96+rax], xmm6
movaps XMMWORD PTR [rdi+112+rax], xmm7
sub rax, -128
cmp rax, 8000000
jne .L2
rep ret
第四次测试
g++ -O3 -fopenmp -mavx test.cpp
./a.out
mul time 0.87 s, 25.3 GB/s, efficiency 74.3%
mul_SSE time 0.88 s, 24.9 GB/s, efficiency 73.0%
mul_SSE_NT time 1.07 s, 20.6 GB/s, efficiency 60.5%
mul_SSE_OMP time 0.76 s, 29.0 GB/s, efficiency 85.2%
现在非内在函数是最快的(不包括 OpenMP 版本)。
因此在这种情况下没有理由使用内在函数或内联汇编,因为我们可以通过适当的编译器选项(例如-O3
、-funroll-loops
、-mavx
)获得最佳性能。
测试系统:Ubuntu 16.10,Skylake (i7-6700HQ@2.60GHz),32GB RAM。最大内存带宽(34.1 GB/s)https://ark.intel.com/products/88967/Intel-Core-i7-6700HQ-Processor-6M-Cache-up-to-3_50-GHz
这是另一个值得考虑的解决方案。 The cmp
instruction is not necessary 如果我们从 -N 数到零并以N+i
访问数组。 GCC 早就应该解决这个问题了。它消除了一条指令(尽管由于宏操作融合, cmp 和 jmp 通常算作一个微操作)。
void mul_SSE_v2(double* a, double* b)
for (ptrdiff_t i = -N; i<0; i+=2)
_mm_store_pd(&a[N + i], _mm_mul_pd(_mm_load_pd(&a[N + i]),_mm_load_pd(&b[N + i])));
与-O3
组装
mul_SSE_v2(double*, double*):
mov rax, -1000000
.L9:
movapd xmm0, XMMWORD PTR [rdi+8000000+rax*8]
mulpd xmm0, XMMWORD PTR [rsi+8000000+rax*8]
movaps XMMWORD PTR [rdi+8000000+rax*8], xmm0
add rax, 2
jne .L9
rep ret
这种优化只可能有助于数组适合,例如L1 缓存,即不从主存读取。
我终于找到了一种方法来获得不生成 cmp
指令的普通 C 解决方案。
void mul_v2(aligned_double* __restrict a, aligned_double* __restrict b)
for (int i = -N; i<0; i++) a[i] *= b[i];
然后从像 mul_v2(&a[N],&b[N])
这样的单独对象文件中调用函数,所以这可能是最好的解决方案。但是,如果您从与在 GCC 中定义的对象文件(翻译单元)相同的对象文件(翻译单元)中调用该函数,则会再次生成 cmp
指令。
还有,
void mul_v3(aligned_double* __restrict a, aligned_double* __restrict b)
for (int i = -N; i<0; i++) a[N+i] *= b[N+i];
仍会生成cmp
指令并生成与mul
函数相同的程序集。
函数mul_SSE_NT
很傻。它使用非临时存储,仅在仅写入内存时才有用,但由于函数读取和写入同一地址,非临时存储不仅无用,而且结果较差。
此答案的先前版本获得了错误的带宽。原因是数组没有初始化。
【讨论】:
我在我的系统上尝试了您的程序以及我的 2 个实现,并将结果添加到开场问题中。我非常喜欢这个答案,因为它非常详细,并且还提供了代码和测量值作为比较,尽管在我接受这个作为答案之前,我只想澄清一下问题本身。为什么普通的 C/C++ 实现运行在 1.33(在我的系统上),而 SIMD 实现运行在 1.09?这是因为它受内存限制吗?如果是,如何知道您的程序何时受内存限制?有什么办法可以优化吗? @fighting_falcon93,因为您的操作受内存带宽限制,因此它不会随 SIMD 通道数或线程数而扩展。但是,它仍然可以从多线程、展开和 SIMD 中受益。这是大多数人不欣赏的部分。我从一开始就更新了我的答案,提供了更多细节。 @fighting_falcon93 我忘了回答你关于 OpenMP 的问题。如果您使用-fopenmp
编译,您将看到call GOMP_parallel
和其他代码,因此OpenMP 程序集与没有godbolt.org/g/yZkH23 的情况不同。
@fighting_falcon93,我修正了我的答案。问题是我使用的是未初始化的数组。 memset(b, 1, sizeof *a * N)
已修复!我重写了代码。它现在只是一个文件,而且更干净。我清理了其余的答案。我现在很满意。
@fighting_falcon93,好的,我通过时间校正再次更新了我的答案。让我知道你的想法。我从这个问题中学到了很多东西。【参考方案2】:
我想为这个问题添加另一个观点。 如果没有内存限制,SIMD 指令会大大提高性能。但是在当前示例中,内存加载和存储操作太多,CPU 计算太少。因此 CPU 可以及时处理传入的数据,而无需使用 SIMD。 如果您使用其他类型的数据(例如 32 位浮点数)或更复杂的算法,内存吞吐量不会限制 CPU 性能,使用 SIMD 会带来更多优势。
【讨论】:
那是我最初的想法:内存带宽受限。但在我的测试中,我仍然看到 N=1000000(2 个双数组,因此 16 MB)的矢量化有显着改进。 好吧,考虑到 OP(最后一个实验)中的循环展开实验,我认为我们可以得出结论,CPU 根本无法并行执行物理上可能的所有内存提取。因此,OP 已经达到了内存障碍,只是不是在吞吐量方面,而是在延迟方面。 @Ermlg 好点。有什么方法可以确定实现是否受内存限制?或任何其他类型的界限,例如分支错误预测界限或输入/输出界限?【参考方案3】:你的 asm 代码真的没问题。 不是你衡量它的方式。 正如我在 cmets 中指出的那样,您应该:
a) 使用更多的迭代次数 - 100 万次对于现代 CPU 来说不算什么
b) 使用 HPT 进行测量
c) 使用 RDTSC 或 RDTSCP 计算实际 CPU 时钟数
另外你为什么害怕-O3 opt?不要忘记为您的平台构建代码,因此请使用 -march=native。如果您的 CPU 支持 AVX 或 AVX2 编译器将借此机会生成更好的代码。
接下来 - 如果你知道你的代码,给编译器一些关于别名和对齐的提示。
这是我的 mul_c
版本 - 是的,它是 GCC 特定的,但你表明你使用了 GCC
void mul_c(double* restrict a, double* restrict b)
a = __builtin_assume_aligned (a, 16);
b = __builtin_assume_aligned (b, 16);
for (int i = 0; i != 1000000; ++i)
a[i] = a[i] * b[i];
它将产生:
mul_c(double*, double*):
xor eax, eax
.L2:
movapd xmm0, XMMWORD PTR [rdi+rax]
mulpd xmm0, XMMWORD PTR [rsi+rax]
movaps XMMWORD PTR [rdi+rax], xmm0
add rax, 16
cmp rax, 8000000
jne .L2
rep ret
如果你有 AVX2 并确保数据是 32 字节对齐的,它将变成
mul_c(double*, double*):
xor eax, eax
.L2:
vmovapd ymm0, YMMWORD PTR [rdi+rax]
vmulpd ymm0, ymm0, YMMWORD PTR [rsi+rax]
vmovapd YMMWORD PTR [rdi+rax], ymm0
add rax, 32
cmp rax, 8000000
jne .L2
vzeroupper
ret
因此,如果编译器可以为您完成,则无需手工制作 asm ;)
【讨论】:
我尝试用 RDTSCP 测量运行时间,但我用新的代码和结果更新了我的问题。正如我在更新中所写,周期的数量变化很大,这可能是因为操作系统在其间运行其他进程。有没有办法只计算我的程序期间的周期?另外,为什么mul_asm2
计算周期的速度这么快?我不使用-O3
的原因是因为我稍后将运行代码的系统不允许我指定标志,它使用-O2
,否则我会使用-O3
;)另外,谢谢你的提示,我不知道这样的提示是可能的。
另外,我稍后运行它的系统支持 AVX2,但我现在正在使用的系统不支持,所以我只使用 128 位 (XMM)现在注册。稍后我会将其更改为 256 位寄存器 (YMM)。将 AVX-512 与 512 位寄存器 (ZMM) 一起使用会很酷,但两个系统都不支持它:'(
@fighting_falcon93 修复 C 源代码而不是编写 asm 的要点是,您可以在不更改源代码的情况下为两个系统编译(在您的系统上,它会在没有 AVX2 的情况下编译,在目标上它将使用AVX2(如果使用了正确的编译时间开关))。那么,如果 C 足以生成最佳矢量化代码,为什么还要修复 asm?
@Ped7g 主要是想学习。我认为编写汇编并击败编译器很有趣,而且我经常注意到编译器会做一些没有完全优化的愚蠢事情。我做了很多编程,其中性能非常重要,每毫秒越少越好,并且您希望代码尽可能快地运行,例如在游戏中以及与其他网站上拥有更快代码的人竞争时CodeChef 等所以我正在尝试寻找新的方法来将我的实现的性能推到极限:)以上是关于为啥这个 SIMD 乘法不比非 SIMD 乘法快?的主要内容,如果未能解决你的问题,请参考以下文章
如何使用内部函数 C++ 将 3 个加法和 1 个乘法转换为矢量化 SIMD
这个 Delphi 6 位图修改代码可以用 SIMD 或其他方法加速吗?