优化的2x2矩阵乘法:慢速装配与快速SIMD

Question

优化的2x2矩阵乘法:慢速装配与快速SIMD

问题

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

到目前为止我看过的是什么

所有相关论文,课程网页,许多SO Q&As处理主题(太多无法列出),我已经在我的计算机上编译了OpenBLAS,查看了OpenBLAS,GotoBLAS和BLIS源代码,Agner的手册.

硬件

我的CPU是Intel i5 - 540M.您可以在cpu-world.com上找到相关的CPUID信息.微体系结构是Nehalem(westmere),因此理论上每循环每个核心可以计算4个双精度触发器.我将只使用一个核心(没有OpenMP),所以对于超线程关闭和4步Intel Turbo Boost,我应该会看到一个高峰( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops.作为参考,两个核心都运行在峰值,英特尔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;
}

Run Code Online (Sandbox Code Playgroud)

我使用朴素矩阵矩阵乘法(如下所示)计算一个真正的答案,它允许我通过视觉或通过计算所有元素的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];

Run Code Online (Sandbox Code Playgroud)

为了运行代码并获得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);

Run Code Online (Sandbox Code Playgroud)

SIMD代码

我的CPU支持128位向量,所以我可以double在每个向量中适应2 s.这是我在内核中进行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);
}

Run Code Online (Sandbox Code Playgroud)

Assmebly(英特尔语法)

我的第一次尝试是为这个部分创建一个单独的汇编例程,并从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]                   \n\t"
    "mov        ecx, %[B]                   \n\t"
    "mov        eax, %[C]                   \n\t"
    "movapd     xmm3, XMMWORD PTR [ecx]     \n\t"
    "movapd     xmm2, XMMWORD PTR [ecx+16]  \n\t"
    "movddup    xmm1, QWORD PTR [edx]       \n\t"
    "mulpd      xmm1, xmm3                  \n\t"
    "addpd      xmm1, XMMWORD PTR [eax]     \n\t"
    "movddup    xmm0, QWORD PTR [edx+8]     \n\t"
    "mulpd      xmm0, xmm2                  \n\t"
    "addpd      xmm0, xmm1                  \n\t"
    "movapd     XMMWORD PTR [eax], xmm0     \n\t"
    "movddup    xmm4, QWORD PTR [edx+16]    \n\t"
    "mulpd      xmm4, xmm3                  \n\t"
    "addpd      xmm4, XMMWORD PTR [eax+16]  \n\t"
    "movddup    xmm5, QWORD PTR [edx+24]    \n\t"
    "mulpd      xmm5, xmm2                  \n\t"
    "addpd      xmm5, xmm4                  \n\t"
    "movapd     XMMWORD PTR [eax+16], xmm5  \n\t"
    : // no outputs 
    : // inputs
    [A] "m" (A),
    [B] "m" (B), 
    [C] "m" (C)
    : //register clobber
    "memory",
    "edx","ecx","eax",
    "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5"
    );
}

Run Code Online (Sandbox Code Playgroud)

结果

我用以下标志编译我的代码:

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

Run Code Online (Sandbox Code Playgroud)

结果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

Run Code Online (Sandbox Code Playgroud)

如果我强制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

Run Code Online (Sandbox Code Playgroud)

问题

如果内联ASM和SIMD实现产生相同的汇编输出,为什么汇编版本要慢得多？就像内联汇编没有内联一样,第二组结果显示"内联"ASM与"无内线"SIMD的性能相同.我能找到的唯一解释是在Agner Fog第2卷第6页:

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

但两个版本的汇编输出完全相同.
为什么我在第一组结果中看到了44个Gflops？这远远高于我计算的12 Gflops峰值,如果我使用单精度计算运行两个核心,这就是我所期望的.

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

Answer 1

Z b*_*son 3

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

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

Run Code Online (Sandbox Code Playgroud)

如果没有这条线，在 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 %f\n", C[0], C[1], C[2], C[3]);
  //gflops = (double) (2.0*n*n*n)/time3/1.0e9*maxiter;
  printf("time %f\n", dtime);
}

Run Code Online (Sandbox Code Playgroud)

归档时间：	11 年，5 月前
查看次数：	1860 次
最近记录：	11 年，5 月前