向量的点积与SIMD

4
我正在尝试使用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代码速度会快8倍左右。你们有什么想法吗?
感谢大家的帮助 :)

1
时间复杂度是多少?而且使用malloc来分配tmp是没有意义的。这样做比直接使用栈要慢得多。 - Sami Kuhmonen
1
你检查过编译器是否自动向量化标量代码了吗? - EOF
相关:如何优化我的AVX点积实现?展示了如何将hsum移出循环。为什么Haswell上的mulss只需要3个周期,与Agner的指令表不同?(使用多个累加器展开FP循环)展示了如何通过多个累加器进一步优化。 - Peter Cordes
AVX2:计算512个浮点数组的点积具有使用内部函数手动向量化的良好循环。 - Peter Cordes
2个回答

14

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

是的,您现在在您“优化”的代码中使用了AVX指令。但是,您还编写了一些编译器现在无法展开的代码,除了普通的向量化之外。

为了比较,让我们看一下编译器实际从您的“慢速”C实现中生成了什么,只是热循环而没有脚注。

ICC,使用-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更为悲观,并将其展开为以下内容:

.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也使用了相同的选项,实际上只能做到和你的“优化”解决方案一样好:

.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 在该循环中还添加了一个头文件,因此它可以使用 vmovaps(对齐内存访问)而不是 vmovups(非对齐内存访问)进行第一次加载操作,更为聪明。


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

ICC:

..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

Clang:

.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

GCC:

.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

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


1
@PaulR 对,是我的错。但这并不会改变编译器优化的选择。两个编译器仍然会进行向量化、展开和交错处理。 - Ext3h
2
在我看来,Clang做得最好。英特尔编译器会产生很多针对非常特定的迭代计数进行优化的代码路径,这样会产生太多的膨胀,不符合我的口味。虽然它确实使用了更少的寄存器,为寄存器重命名和因此的推测执行留下了更多的余地,但我更喜欢Clang的输出,因为它在掩盖指令延迟方面做得更好,尤其是在非英特尔处理器上。 - Ext3h
1
@Ext3h:并不是说-march=ivybridge缺少未对齐的256b加载,而是gcc决定不使用它,因为如果数据在运行时实际上未对齐(而不仅仅是在编译时不知道始终对齐),那么使用vmovups / vinsertf128会提高性能。tune = generic和tune = sandybridge / ivybridge(以及一些AMD调整)有意启用-mavx256-split-unaligned-load这对于-mtune = generic来说可能不是理想的 - Peter Cordes
2
拥有更多的操作可以帮助掩盖有限的内存吞吐量。不,多个累加器隐藏了vfmaddps的延迟,直到您在负载吞吐量(每个时钟最多2个YMM负载)上瓶颈而不是循环传递FMA延迟(Haswell上为5个周期)。使用较少的寄存器,留下更多的余地进行寄存器重命名,从而进行推测执行:重写相同的体系结构寄存器仍需要一个新的PRF条目。保留一些体系结构寄存器对乱序窗口大小几乎没有影响。 - Peter Cordes
2
@PeterCordes 谢谢,已添加链接并替换了有关额外累加器的虚假解释。 - Ext3h
显示剩余8条评论

2

不幸的是,点积算法是一个内存绑定算法(计算次数少于所需内存吞吐量)。因此,即使使用AVX(或AVX2),您也无法有效地执行它。我还以类似的方式实现了这个算法,但我只达到了性能提高60%的水平。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接