我在下面整理了一些示例代码,以说明您如何看待 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
查看生成的代码很有趣:
Scalar:
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(2x 展开):
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 循环每次迭代当然处理四个元素,而标量循环每次迭代仅处理一个元素。对于第三个展开循环,我们有更多指令,但每次迭代处理 8 个元素 - 请注意,相对于没有循环展开的 SIMD 循环,循环管理指令的比例已减少。
计时数据是使用 2.6 GHz Core i7 Haswell CPU 在 Mac OS X 10.10 上使用 gcc 4.8 收集的。然而,在任何当前合理的 x86 CPU 上,性能结果应该相似。