c - 带有 SIMD 的 vector 的点积

Question

我正在尝试使用 SIMD 指令来加速我的 C 代码中的点积计算。但是，我的函数的运行时间大致相等。如果有人能解释为什么以及如何加速计算，那就太好了。

具体来说，我试图计算两个数组的点积，其中大约有 10,000 个元素。我的常规 C 函数如下:

 float my_dotProd( float const * const x, float const * const y, size_t const N ){
   // N is the number of elements in the arrays
   size_t i;
   float out=0;

   for( i=0; i < N; ++i ){
     out += x[i] * y[i];
   }

   return out;
 }

我使用 AVX SIMD 命令的功能如下:

 void my_malloc( size_t nBytes, void ** ptrPtr ){
   int boundary = 32;
   posix_memalign( ptrPtr, boundary, nBytes );
 }

 float cimpl_sum_m128( __m128 x ){
   float out;
   __m128 sum = x;
   sum = _mm_hadd_ps( sum, sum );
   sum = _mm_hadd_ps( sum, sum );
   out = _mm_cvtss_f32( sum );
   return out;
 }

 float my_sum_m256( __m256 x ){
   float out1, out2;
   __m128 hi = _mm256_extractf128_ps(x, 1);
   __m128 lo = _mm256_extractf128_ps(x, 0);
   out1 = cimpl_sum_m128( hi );
   out2 = cimpl_sum_m128( lo );
   return out1 + out2;
 }

 float my_dotProd( float const * const x, float const * const y, size_t const N ){
   // N is the number of elements in the arrays
   size_t i=0;
   float out=0;
   float *tmp;

   __m256 summed, *l, *r;

   if( N > 7 ){
     my_malloc( sizeof(float) * 8, (void**) &tmp );
     summed = _mm256_set1_ps(0.0f);
     l = (__m256*) x;
     r = (__m256*) y;

     for( i=0; i < N-7; i+=8, ++l, ++r ){
       summed = _mm256_add_ps( summed, _mm256_mul_ps( *l, *r ) );
     }
     _mm256_store_ps( tmp, summed );

     out += my_sum_m256( summed );
     free( tmp );
   }

   for( ; i < N; ++i ){
     out += x[i] * y[i];
   }

   return out;
 }

我的测试程序是:

 int test_dotProd(){
   float *x, *y;
   size_t i, N;
   float answer, result;
   float err;

   N = 100000;  // Fails

   my_malloc( sizeof(float) * N, (void**) &x );
   my_malloc( sizeof(float) * N, (void**) &y );

   answer = 0;
   for( i=0; i<N; ++i ){
     x[i]=i; y[i]=i;
     answer += (float)i * (float)i;
   }

   result = my_dotProd( x, y, N );

   err = fabs( result - answer ) / answer;

   free( x );
   free( y );
   return err < 5e-7;
 }

我正在使用时钟来测量这样的运行时间:

 timeStart = clock();
 testStatus = test_dotProd();
 timeTaken = (int)( clock() - timeStart );

我意识到 my_sum_m256 操作可以提高效率，但我认为这对运行时的影响应该很小。我猜想 SIMD 代码的速度会快八倍。有什么想法吗？

谢谢大家的帮助 :)

Answer 1

首先:您不应该假设您可以比编译器更好地优化。

是的，您现在正在“优化”代码中使用 AVX 指令。但是您还编写了编译器现在无法展开的代码，除了简单的矢量化之外。

为了进行比较，让我们看看编译器实际上会从您的“慢”C 实现中生成什么，只是没有页脚的热循环。

ICC, compiled with -O3 -march=skylake -ffast-math :

..B1.13:
    vmovups   ymm2, YMMWORD PTR [rsi+rdi*4]
    vmovups   ymm3, YMMWORD PTR [32+rsi+rdi*4]
    vfmadd231ps ymm1, ymm2, YMMWORD PTR [r8+rdi*4]
    vfmadd231ps ymm0, ymm3, YMMWORD PTR [32+r8+rdi*4]
    add       rdi, 16
    cmp       rdi, rax
    jb        ..B1.13

Clang, with the same parameters更加悲观，并将其展开为以下内容:

.LBB0_4:
    vmovups ymm4, ymmword ptr [rsi + 4*rcx]
    vmovups ymm5, ymmword ptr [rsi + 4*rcx + 32]
    vmovups ymm6, ymmword ptr [rsi + 4*rcx + 64]
    vmovups ymm7, ymmword ptr [rsi + 4*rcx + 96]
    vfmadd132ps     ymm4, ymm0, ymmword ptr [rdi + 4*rcx]
    vfmadd132ps     ymm5, ymm1, ymmword ptr [rdi + 4*rcx + 32]
    vfmadd132ps     ymm6, ymm2, ymmword ptr [rdi + 4*rcx + 64]
    vfmadd132ps     ymm7, ymm3, ymmword ptr [rdi + 4*rcx + 96]
    vmovups ymm0, ymmword ptr [rsi + 4*rcx + 128]
    vmovups ymm1, ymmword ptr [rsi + 4*rcx + 160]
    vmovups ymm2, ymmword ptr [rsi + 4*rcx + 192]
    vmovups ymm3, ymmword ptr [rsi + 4*rcx + 224]
    vfmadd132ps     ymm0, ymm4, ymmword ptr [rdi + 4*rcx + 128]
    vfmadd132ps     ymm1, ymm5, ymmword ptr [rdi + 4*rcx + 160]
    vfmadd132ps     ymm2, ymm6, ymmword ptr [rdi + 4*rcx + 192]
    vfmadd132ps     ymm3, ymm7, ymmword ptr [rdi + 4*rcx + 224]
    add     rcx, 64
    add     rax, 2
    jne     .LBB0_4

令人惊讶的是，两个编译器都已经能够使用 AVX 指令，不需要内在的黑客攻击。

但更有趣的是，两个编译器都认为一个累加寄存器不足以使 AVX 管道饱和，而是分别使用了 2 个和 4 个累加寄存器。进行更多操作有助于屏蔽 FMA 的延迟，直至达到实际内存吞吐量限制。

只是不要忘记 -ffast-math编译器选项，如果没有从矢量化循环中提取最终累加是不合法的。

GCC, also with the same options , 实际上“仅”与您的“优化”解决方案一样好:

.L7:
    add     r8, 1
    vmovaps ymm3, YMMWORD PTR [r9+rax]
    vfmadd231ps     ymm1, ymm3, YMMWORD PTR [rcx+rax]
    add     rax, 32
    cmp     r8, r10
    jb      .L7

然而，GCC 在向该循环添加 header 方面仍然更聪明一点，因此它可以使用 vmovaps (对齐的内存访问)而不是 vmovups (未对齐的内存访问)用于第一次加载。

为了完整起见，使用纯 AVX ( -O3 -march=ivybridge -ffast-math ):

国际刑事法院:

..B1.12:
    vmovups   xmm2, XMMWORD PTR [r8+rdi*4]
    vmovups   xmm5, XMMWORD PTR [32+r8+rdi*4]
    vinsertf128 ymm3, ymm2, XMMWORD PTR [16+r8+rdi*4], 1
    vinsertf128 ymm6, ymm5, XMMWORD PTR [48+r8+rdi*4], 1
    vmulps    ymm4, ymm3, YMMWORD PTR [rsi+rdi*4]
    vmulps    ymm7, ymm6, YMMWORD PTR [32+rsi+rdi*4]
    vaddps    ymm1, ymm1, ymm4
    vaddps    ymm0, ymm0, ymm7
    add       rdi, 16
    cmp       rdi, rax
    jb        ..B1.12

铛:

.LBB0_5:
    vmovups xmm4, xmmword ptr [rdi + 4*rcx]
    vmovups xmm5, xmmword ptr [rdi + 4*rcx + 32]
    vmovups xmm6, xmmword ptr [rdi + 4*rcx + 64]
    vmovups xmm7, xmmword ptr [rdi + 4*rcx + 96]
    vinsertf128     ymm4, ymm4, xmmword ptr [rdi + 4*rcx + 16], 1
    vinsertf128     ymm5, ymm5, xmmword ptr [rdi + 4*rcx + 48], 1
    vinsertf128     ymm6, ymm6, xmmword ptr [rdi + 4*rcx + 80], 1
    vinsertf128     ymm7, ymm7, xmmword ptr [rdi + 4*rcx + 112], 1
    vmovups xmm8, xmmword ptr [rsi + 4*rcx]
    vmovups xmm9, xmmword ptr [rsi + 4*rcx + 32]
    vmovups xmm10, xmmword ptr [rsi + 4*rcx + 64]
    vmovups xmm11, xmmword ptr [rsi + 4*rcx + 96]
    vinsertf128     ymm8, ymm8, xmmword ptr [rsi + 4*rcx + 16], 1
    vmulps  ymm4, ymm8, ymm4
    vaddps  ymm0, ymm4, ymm0
    vinsertf128     ymm4, ymm9, xmmword ptr [rsi + 4*rcx + 48], 1
    vmulps  ymm4, ymm4, ymm5
    vaddps  ymm1, ymm4, ymm1
    vinsertf128     ymm4, ymm10, xmmword ptr [rsi + 4*rcx + 80], 1
    vmulps  ymm4, ymm4, ymm6
    vaddps  ymm2, ymm4, ymm2
    vinsertf128     ymm4, ymm11, xmmword ptr [rsi + 4*rcx + 112], 1
    vmulps  ymm4, ymm4, ymm7
    vaddps  ymm3, ymm4, ymm3
    add     rcx, 32
    cmp     rax, rcx
    jne     .LBB0_5

海湾合作委员会:

.L5:
    vmovups xmm3, XMMWORD PTR [rdi+rax]
    vinsertf128     ymm1, ymm3, XMMWORD PTR [rdi+16+rax], 0x1
    vmovups xmm4, XMMWORD PTR [rsi+rax]
    vinsertf128     ymm2, ymm4, XMMWORD PTR [rsi+16+rax], 0x1
    add     rax, 32
    vmulps  ymm1, ymm1, ymm2
    vaddps  ymm0, ymm0, ymm1
    cmp     rax, rcx
    jne     .L5

应用了几乎相同的优化，只是一些额外的操作，因为缺少 FMA 并且不建议对 Ivy Bridge 使用未对齐的 256 位加载。