优化的2x2矩阵乘法：与快速SIMD相比，组装速度慢

2018-06-06 16:02:26

问题

我正在研究高性能的矩阵乘法算法，如OpenBLAS或GotoBLAS，我正试图重现一些结果。这个问题涉及矩阵乘法算法的内核。具体来说，我在计算C += AB ，其中A和B是我的CPU的峰值速度下的double类型的2×2矩阵。有两种方法可以做到这一点。一种方法是使用SIMD指令。第二种方法是使用SIMD寄存器在汇编中直接编码。

迄今为止我所看到的

所有相关论文，当然网页，我的电脑上很多很多的SO Q＆作为处理对象（举不胜举），我已经编译OpenBLAS，通过OpenBLAS，GotoBLAS，以及BLIS源代码，昂纳手册看了。

硬件

我的CPU是Intel i5-540M。您可以在cpu-world.com上找到相关的CPUID信息。微体系结构是Nehalem（westmere），所以它在理论上可以计算每个内核每个周期4个双精度触发器。我将仅使用一个芯（无OpenMP的），所以与关闭超线程和4步英特尔智能加速，我应该看到的峰值( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops 。作为参考，两个内核运行在最高峰时，Intel Turbo Boost提供了两步加速，我应该得到22.4 DP Gflops的理论峰值。

建立

我将我的2x2矩阵声明为double ，并用随机条目初始化它们，如下面的代码片段所示。

srand(time(NULL));
const int n = 2;
double A[n*n];
double B[n*n];
double C[n*n];
double T[n*n];
for(int i = 0; i < n*n; i++){
    A[i] = (double) rand()/RAND_MAX;
    B[i] = (double) rand()/RAND_MAX;
    C[i] = 0.0;
}

我使用天真的矩阵 - 矩阵乘法（下面显示）来计算一个真正的答案，这允许我通过视觉或计算所有元素的L2范数来检查我的结果

// "true" answer
for(int i = 0; i < n; i++)
    for(int j = 0; j < n; j++)
        for(int k = 0; k < n; k++)
            T[i*n + j] += A[i*n + k]*B[k*n + j];

运行代码并获得GFLOPS的估计，我呼吁每个乘法函数一次热身，然后执行它里面for循环的maxiter次，确保到零C ，因为我计算每次矩阵C += AB 。 for循环放置在两个clock()语句中，并用于估计Gflops。代码片段打击说明了这一部分。

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0;
mult2by2(A,B,C); //warmup
time1 = clock();
for(int i = 0; i < maxiter; i++){
        mult2by2(A,B,C);
        C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0;
}
time2 = clock() - time1;
time3 = (double)(time2)/CLOCKS_PER_SEC;
gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter;
mult2by2(A,B,C); // to compute the norm against T
norm = L2norm(n,C,T);

SIMD代码

我的CPU支持128位向量，所以我可以在每个向量中放置2个double精度值。这是我在内核中进行2x2矩阵乘法的主要原因。 SIMD代码一次计算整行C

    inline void 
    __attribute__ ((gnu_inline))        
    __attribute__ ((aligned(16))) mult2by2B(        
            const double* restrict A,
            const double* restrict B,
            double* restrict C
        )

    {

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4;
    xmm0 = _mm_load_pd(C);
    xmm1 = _mm_load1_pd(A);
    xmm2 = _mm_load_pd(B);
    xmm3 = _mm_load1_pd(A + 1);
    xmm4 = _mm_load_pd(B + 2);
    xmm1 = _mm_mul_pd(xmm1,xmm2);
    xmm2 = _mm_add_pd(xmm1,xmm0);
    xmm1 = _mm_mul_pd(xmm3,xmm4);
    xmm2 = _mm_add_pd(xmm1,xmm2);
    _mm_store_pd(C,xmm2);

    xmm0 = _mm_load_pd(C + 2);
    xmm1 = _mm_load1_pd(A + 2);
    xmm2 = _mm_load_pd(B);
    xmm3 = _mm_load1_pd(A + 3);
    //xmm4 = _mm_load_pd(B + 2);
    xmm1 = _mm_mul_pd(xmm1,xmm2);
    xmm2 = _mm_add_pd(xmm1,xmm0);
    xmm1 = _mm_mul_pd(xmm3,xmm4);
    xmm2 = _mm_add_pd(xmm1,xmm2);
    _mm_store_pd(C + 2,xmm2);
}

（英特尔语法）

我的第一个尝试是为这部分创建一个单独的汇编例程，并从main例程中调用它。然而，它非常缓慢，因为我不能内联extern函数。我将程序集编写为内联程序集，如下所示。它是相同的，这是由生产gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel 。根据我所了解的Nehalem微架构图，该处理器可以并行执行SSE ADD ， SSE MUL和SSE MOV ，这解释了MUL ， ADD ， MOV指令的交错。您会注意到上述SIMD指令的顺序不同，因为我对Agner Fog的手册有不同的理解。尽管如此， gcc很聪明，上面的SIMD代码编译为内联版本中显示的程序集。

inline void 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(16))) mult2by2A
    (   
        const double* restrict A,
        const double* restrict B,
        double* restrict C
    )
    {
    __asm__ __volatile__
    (
    "mov        edx, %[A]                   nt"
    "mov        ecx, %[B]                   nt"
    "mov        eax, %[C]                   nt"
    "movapd     xmm3, XMMWORD PTR [ecx]     nt"
    "movapd     xmm2, XMMWORD PTR [ecx+16]  nt"
    "movddup    xmm1, QWORD PTR [edx]       nt"
    "mulpd      xmm1, xmm3                  nt"
    "addpd      xmm1, XMMWORD PTR [eax]     nt"
    "movddup    xmm0, QWORD PTR [edx+8]     nt"
    "mulpd      xmm0, xmm2                  nt"
    "addpd      xmm0, xmm1                  nt"
    "movapd     XMMWORD PTR [eax], xmm0     nt"
    "movddup    xmm4, QWORD PTR [edx+16]    nt"
    "mulpd      xmm4, xmm3                  nt"
    "addpd      xmm4, XMMWORD PTR [eax+16]  nt"
    "movddup    xmm5, QWORD PTR [edx+24]    nt"
    "mulpd      xmm5, xmm2                  nt"
    "addpd      xmm5, xmm4                  nt"
    "movapd     XMMWORD PTR [eax+16], xmm5  nt"
    : // no outputs 
    : // inputs
    [A] "m" (A),
    [B] "m" (B), 
    [C] "m" (C)
    : //register clobber
    "memory",
    "edx","ecx","eax",
    "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5"
    );
}

结果

我用下面的标志编译我的代码：

gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel

maxiter = 1000000000的结果如下：

********** Inline ASM
L2 norm: 0.000000e+000, Avg. CPU time: 9.563000, Avg. Gflops: 1.673115

********** SIMD Version
L2 norm: 0.000000e+000, Avg. CPU time: 0.359000, Avg. Gflops: 44.568245

如果我强制SIMD版本不用__attribute__ ((noinline))内联，结果是：

********** Inline ASM
L2 norm: 0.000000e+000, Avg. CPU time: 11.155000, Avg. Gflops: 1.434334

********** SIMD Version
L2 norm: 0.000000e+000, Avg. CPU time: 11.264000, Avg. Gflops: 1.420455

问题

如果内联ASM和SIMD实现产生相同的汇编输出，为什么汇编版本要慢得多？就好像内联程序集没有被内联，而第二组结果显示“内联”ASM与“非内联”SIMD的性能完全相同。我能找到的唯一解释是在Agner Fog第2卷第6页：

编译代码可能比汇编代码更快，因为编译器可以进行程序间优化和整个程序优化。汇编程序员通常必须使用定义良好的调用接口来定义明确的函数，该接口遵从所有调用约定，以使代码可测试和可验证。这可以防止许多编译器使用的优化方法，例如函数内联，寄存器分配，常量传播，跨函数的通用子表达式消除，跨函数调度等。这些优点可以通过使用具有内部函数而不是汇编代码的C ++代码来获得。

但两个版本的汇编器输出完全相同。

为什么我在第一组结果中看到44 Gflops？这比我计算出来的12 Gflops峰值高出一些，而且如果我用单精度计算运行两个内核，这是我所期望的。

编辑1评论说可能有死代码消除我可以确认这是发生在SIMd指令。 -S输出显示SIMD唯一零C矩阵的for循环。我可以通过关闭-O0编译器优化来禁用它。在这种情况下，SIMD的运行速度是ASM的3倍，但ASM仍以完全相同的速度运行。现在的规范也不为零，但在10 ^ -16时仍然可以。我还看到内联ASM版本正在使用APP和NO_APP标签进行内联，但它也在for循环中展开了8次。我认为多次展开会严重影响性能，因为我通常展开循环4次。根据我的经验，任何事情似乎都会降低性能。

GCC正在使用内部函数mult2by2B优化您的内联函数，因为该行

C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0;

如果没有这条线，Coliru的电脑需要2.9秒http://coliru.stacked-crooked.com/a/992304f5f672e257

与该行，它只需要0.000001 http://coliru.stacked-crooked.com/a/9722c39bb6b8590a

你也可以在程序集中看到它。如果你把下面的代码放到http://gcc.godbolt.org/中，你会看到用那行代码完全跳过了这个函数。

但是，当您将内联程序集GCC不优化函数mult2by2A ，即使它将内联函数内联。你也可以在程序集中看到它。

#include <stdio.h>
#include <emmintrin.h>                 // SSE2
#include <omp.h>

inline void 
    __attribute__ ((gnu_inline))        
    __attribute__ ((aligned(16))) mult2by2B(        
            const double* __restrict A,
            const double* __restrict B,
            double* __restrict C
        )

    {

    register __m128d xmm0, xmm1, xmm2, xmm3, xmm4;
    xmm0 = _mm_load_pd(C);
    xmm1 = _mm_load1_pd(A);
    xmm2 = _mm_load_pd(B);
    xmm3 = _mm_load1_pd(A + 1);
    xmm4 = _mm_load_pd(B + 2);
    xmm1 = _mm_mul_pd(xmm1,xmm2);
    xmm2 = _mm_add_pd(xmm1,xmm0);
    xmm1 = _mm_mul_pd(xmm3,xmm4);
    xmm2 = _mm_add_pd(xmm1,xmm2);
    _mm_store_pd(C,xmm2);

    xmm0 = _mm_load_pd(C + 2);
    xmm1 = _mm_load1_pd(A + 2);
    xmm2 = _mm_load_pd(B);
    xmm3 = _mm_load1_pd(A + 3);
    //xmm4 = _mm_load_pd(B + 2);
    xmm1 = _mm_mul_pd(xmm1,xmm2);
    xmm2 = _mm_add_pd(xmm1,xmm0);
    xmm1 = _mm_mul_pd(xmm3,xmm4);
    xmm2 = _mm_add_pd(xmm1,xmm2);
    _mm_store_pd(C + 2,xmm2);
}

int main() {
  double A[4], B[4], C[4];
  int maxiter = 10000000;
  //int maxiter = 1000000000;
  double dtime;
  dtime = omp_get_wtime();
  for(int i = 0; i < maxiter; i++){
        mult2by2B(A,B,C);
        C[0] = 0.0; C[1] = 0.0; C[2] = 0.0; C[3] = 0.0;
  }
  dtime = omp_get_wtime() - dtime;
  printf("%f %f %f %fn", C[0], C[1], C[2], C[3]);
  //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter;
  printf("time %fn", dtime);
}

链接地址: http://www.djcxy.com/p/20633.html

上一篇: Optimized 2x2 matrix multiplication: Slow assembly versus fast SIMD

下一篇: How to download Google Chart API?