我试图了解使用 SIMD 矢量化的好处，并编写了一个简单的演示程序代码来查看利用矢量化 (SIMD) 的算法相对于另一个算法的速度增益是多少。以下是 2 种算法:
Alg_A - 不支持 vector :
#include <stdio.h>
#define SIZE 1000000000
int main() {
printf("Algorithm with NO vector support\n");
int a[] = {1, 2, 3, 4};
int b[] = {5, 6, 7, 8};
int i = 0;
printf("Running loop %d times\n", SIZE);
for (; i < SIZE; i++) {
a[0] = a[0] + b[0];
a[1] = a[1] + b[1];
a[2] = a[2] + b[2];
a[3] = a[3] + b[3];
}
printf("A: [%d %d %d %d]\n", a[0], a[1], a[2], a[3]);
}
Alg_B - 支持 vector :
#include <stdio.h>
#define SIZE 1000000000
typedef int v4_i __attribute__ ((vector_size(16)));
union Vec4 {
v4_i v;
int i[4];
};
int main() {
printf("Algorithm with vector support\n\n");
union Vec4 a, b;
a.i[0] = 1, a.i[1] = 2, a.i[2] = 3, a.i[3] = 4;
b.i[0] = 5, b.i[1] = 5, b.i[2] = 7, b.i[3] = 8;
int i = 0;
printf("Running loop %d times\n", SIZE);
for (; i < SIZE; i++) {
a.v = a.v + b.v;
}
printf("A: [%d %d %d %d]\n", a.i[0], a.i[1], a.i[2], a.i[3]);
}
编译完成如下:
算法A:
gcc -ggdb -mno-sse -mno-sse2 -mno-sse3 -mno-sse4 -mno-sse4.1 -mno-sse4.2 -c non_vector_support.c
gcc non_vector_support.o -o non_vector_support
算法B:
gcc -ggdb -c vector_support.c
gcc vector_support.o -o vector_support
查看两种算法的生成代码，我可以看到编译没有做任何像“自动矢量化”这样的技巧(例如，看不到 xmm 寄存器):
算法A:
for (; i < SIZE; i++) {
74: eb 30 jmp a6 <main+0xa6>
a[0] = a[0] + b[0];
76: 8b 55 d8 mov -0x28(%rbp),%edx
79: 8b 45 e8 mov -0x18(%rbp),%eax
7c: 01 d0 add %edx,%eax
7e: 89 45 d8 mov %eax,-0x28(%rbp)
a[1] = a[1] + b[1];
81: 8b 55 dc mov -0x24(%rbp),%edx
84: 8b 45 ec mov -0x14(%rbp),%eax
87: 01 d0 add %edx,%eax
89: 89 45 dc mov %eax,-0x24(%rbp)
a[2] = a[2] + b[2];
8c: 8b 55 e0 mov -0x20(%rbp),%edx
8f: 8b 45 f0 mov -0x10(%rbp),%eax
92: 01 d0 add %edx,%eax
94: 89 45 e0 mov %eax,-0x20(%rbp)
a[3] = a[3] + b[3];
97: 8b 55 e4 mov -0x1c(%rbp),%edx
9a: 8b 45 f4 mov -0xc(%rbp),%eax
9d: 01 d0 add %edx,%eax
9f: 89 45 e4 mov %eax,-0x1c(%rbp)
int a[] = {1, 2, 3, 4};
int b[] = {5, 6, 7, 8};
int i = 0;
printf("Running loop %d times\n", SIZE);
for (; i < SIZE; i++) {
a2: 83 45 d4 01 addl $0x1,-0x2c(%rbp)
a6: 81 7d d4 ff c9 9a 3b cmpl $0x3b9ac9ff,-0x2c(%rbp)
ad: 7e c7 jle 76 <main+0x76>
a[1] = a[1] + b[1];
a[2] = a[2] + b[2];
a[3] = a[3] + b[3];
}
printf("A: [%d %d %d %d]\n", a[0], a[1], a[2], a[3]);
算法B:
for (; i < SIZE; i++) {
74: eb 16 jmp 8c <main+0x8c>
a.v = a.v + b.v;
76: 66 0f 6f 4d d0 movdqa -0x30(%rbp),%xmm1
7b: 66 0f 6f 45 e0 movdqa -0x20(%rbp),%xmm0
80: 66 0f fe c1 paddd %xmm1,%xmm0
84: 0f 29 45 d0 movaps %xmm0,-0x30(%rbp)
union Vec4 a, b;
a.i[0] = 1, a.i[1] = 2, a.i[2] = 3, a.i[3] = 4;
b.i[0] = 5, b.i[1] = 5, b.i[2] = 7, b.i[3] = 8;
int i = 0;
printf("Running loop %d times\n", SIZE);
for (; i < SIZE; i++) {
88: 83 45 cc 01 addl $0x1,-0x34(%rbp)
8c: 81 7d cc ff c9 9a 3b cmpl $0x3b9ac9ff,-0x34(%rbp)
93: 7e e1 jle 76 <main+0x76>
a.v = a.v + b.v;
}
printf("A: [%d %d %d %d]\n", a.i[0], a.i[1], a.i[2], a.i[3]);
当我运行这些程序时，我希望看到速度提高 4 倍，但是对于这种数据大小，增益似乎平均为 =~ 1s，如果将 SIZE 增加到 8000000000 左右，增益是=~ 5s。这是上面代码中值的时间:
算法A:
Algorithm with NO vector support
Running loop 1000000000 times
A: [705032705 1705032706 -1589934589 -589934588]
real 0m3.279s
user 0m3.282s
sys 0m0.000s
算法B:
支持 vector 的算法
Running loop 1000000000 times
A: [705032705 705032706 -1589934589 -589934588]
real 0m2.609s
user 0m2.607s
sys 0m0.004s
计算与循环相关的开销。我针对给定的 SIZE 运行了一个空循环，并获得了 =~ 2.2s 的平均值。这使我的平均速度提高了大约 2.5 倍。
我是否遗漏了代码或编译行中的某些内容？或者，这是否应该是正确的，在这种情况下，如果我在每次迭代中利用 4 个数据点和 1 条指令，为什么有人会知道为什么速度没有提高 4 倍？
提前致谢。
最佳答案
我在下面整理了一些示例代码，以说明您如何看待 SIMD 与标量代码的优势。示例代码有点做作，但要注意的要点是循环中需要有足够的算术运算以减轻加载/存储延迟和循环开销 - 与初始实验一样，单个添加操作是不够的.
此示例将 32 位 int 数据的吞吐量提高了大约 4 倍。 SIMD 循环有两种版本:一种是没有展开的简单循环，另一种是有 2 次展开的替代循环。正如预期的那样，展开的循环要快一些。
#include <assert.h>
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <sys/time.h> // gettimeofday
#include <smmintrin.h> // SSE 4.1
static void foo_scalar(uint32_t *a, const uint32_t *b, const uint32_t *c, size_t n)
{
for (size_t i = 0; i < n; ++i)
{
a[i] = (b[i] + c[i] + 1) / 2;
}
}
static void foo_simd(uint32_t *a, const uint32_t *b, const uint32_t *c, size_t n)
{
size_t i;
#ifndef UNROLL
for (i = 0; i <= n - 4; i += 4)
{
__m128i vb = _mm_loadu_si128((__m128i *)&b[i]);
__m128i vc = _mm_loadu_si128((__m128i *)&c[i]);
__m128i v = _mm_add_epi32(vb, vc);
v = _mm_add_epi32(v, _mm_set1_epi32(1));
v = _mm_srli_epi32(v, 1);
_mm_storeu_si128((__m128i *)&a[i], v);
}
#else
for (i = 0; i <= n - 8; i += 8)
{
__m128i vb0 = _mm_loadu_si128((__m128i *)&b[i]);
__m128i vb1 = _mm_loadu_si128((__m128i *)&b[i + 4]);
__m128i vc0 = _mm_loadu_si128((__m128i *)&c[i]);
__m128i vc1 = _mm_loadu_si128((__m128i *)&c[i + 4]);
__m128i v0 = _mm_add_epi32(vb0, vc0);
__m128i v1 = _mm_add_epi32(vb1, vc1);
v0 = _mm_add_epi32(v0, _mm_set1_epi32(1));
v1 = _mm_add_epi32(v1, _mm_set1_epi32(1));
v0 = _mm_srli_epi32(v0, 1);
v1 = _mm_srli_epi32(v1, 1);
_mm_storeu_si128((__m128i *)&a[i], v0);
_mm_storeu_si128((__m128i *)&a[i + 4], v1);
}
#endif
foo_scalar(&a[i], &b[i], &c[i], n - i);
}
int main(int argc, char *argv[])
{
const size_t kLoops = 100000;
size_t n = 2 * 1024;
struct timeval t0, t1;
double t_scalar_ms, t_simd_ms;
if (argc > 1)
{
n = atoi(argv[1]);
}
printf("kLoops = %zu, n = %zu\n", kLoops, n);
uint32_t * a_scalar = malloc(n * sizeof(uint32_t));
uint32_t * a_simd = malloc(n * sizeof(uint32_t));
uint32_t * b = malloc(n * sizeof(uint32_t));
uint32_t * c = malloc(n * sizeof(uint32_t));
for (size_t i = 0; i < n; ++i)
{
a_scalar[i] = a_simd[i] = 0;
b[i] = rand();
c[i] = rand();
}
gettimeofday(&t0, NULL);
for (size_t k = 0; k < kLoops; ++k)
{
foo_scalar(a_scalar, b, c, n);
}
gettimeofday(&t1, NULL);
t_scalar_ms = ((double)(t1.tv_sec - t0.tv_sec) + (double)(t1.tv_usec - t0.tv_usec) * 1.0e-6) * 1.0e3;
gettimeofday(&t0, NULL);
for (size_t k = 0; k < kLoops; ++k)
{
foo_simd(a_simd, b, c, n);
}
gettimeofday(&t1, NULL);
t_simd_ms = ((double)(t1.tv_sec - t0.tv_sec) + (double)(t1.tv_usec - t0.tv_usec) * 1.0e-6) * 1.0e3;
int64_t sum_scalar = 0, sum_simd = 0;
for (size_t i = 0; i < n; ++i)
{
sum_scalar += a_scalar[i];
sum_simd += a_simd[i];
}
assert(sum_scalar == sum_simd);
printf("t_scalar = %8g ms = %8g ns / point\n", t_scalar_ms, t_scalar_ms / (kLoops * n) * 1e6);
printf("t_simd = %8g ms = %8g ns / point\n", t_simd_ms, t_simd_ms / (kLoops * n) * 1e6);
printf("Speed-up = %2.1fx\n", t_scalar_ms / t_simd_ms);
return 0;
}
编译并运行(没有 SIMD 循环展开):
$ gcc-4.8 -fno-tree-vectorize -std=gnu99 -Wall gros_lalo.c -O3 -msse4.1 && ./a.out
kLoops = 100000, n = 2048
t_scalar = 122.668 ms = 0.598965 ns / point
t_simd = 33.785 ms = 0.164966 ns / point
Speed-up = 3.6x
编译并运行(2x SIMD 循环展开):
$ gcc-4.8 -fno-tree-vectorize -std=gnu99 -Wall gros_lalo.c -O3 -msse4.1 -DUNROLL && ./a.out
kLoops = 100000, n = 2048
t_scalar = 121.897 ms = 0.5952 ns / point
t_simd = 29.07 ms = 0.141943 ns / point
Speed-up = 4.2x
看看生成的代码很有意思:
标量:
xorl %ecx, %ecx
.align 4
L10:
movl 0(%rbp,%rcx,4), %esi
addl (%rbx,%rcx,4), %esi
addl $1, %esi
shrl %esi
movl %esi, (%r15,%rcx,4)
addq $1, %rcx
cmpq %r12, %rcx
jne L10
SIMD(不展开):
xorl %ecx, %ecx
xorl %eax, %eax
.align 4
L18:
movdqu 0(%rbp,%rcx), %xmm2
addq $4, %rax
movdqu (%rbx,%rcx), %xmm1
paddd %xmm2, %xmm1
paddd %xmm3, %xmm1
psrld $1, %xmm1
movdqu %xmm1, (%r14,%rcx)
addq $16, %rcx
cmpq %r9, %rax
jbe L18
SIMD(2 倍展开):
xorl %edx, %edx
xorl %ecx, %ecx
.align 4
L18:
movdqu 0(%rbp,%rdx), %xmm5
addq $8, %rcx
movdqu (%r11,%rdx), %xmm4
movdqu (%rbx,%rdx), %xmm2
movdqu (%r10,%rdx), %xmm1
paddd %xmm5, %xmm2
paddd %xmm4, %xmm1
paddd %xmm3, %xmm2
paddd %xmm3, %xmm1
psrld $1, %xmm2
psrld $1, %xmm1
movdqu %xmm2, 0(%r13,%rdx)
movdqu %xmm1, (%rax,%rdx)
addq $32, %rdx
cmpq %r15, %rcx
jbe L18
请注意，前两个循环中的指令数量相似，但 SIMD 循环当然每次迭代处理四个元素，而标量循环每次迭代仅处理一个元素。对于第三个展开的循环，我们有更多的指令，但我们每次迭代处理八个元素 - 请注意，相对于没有循环展开的 SIMD 循环，循环管理指令的比例已经减少。
时序数据是使用 2.6 GHz Core i7 Haswell CPU 在 Mac OS X 10.10 上使用 gcc 4.8 收集的。然而，性能结果在任何当前合理的 x86 CPU 上应该是相似的。
关于c - 在禁用优化的情况下，演示代码无法显示 4 倍快的 SIMD 速度，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/27075540/