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

Problem

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

What I have looked at so far

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

Hardware

我的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 的理论峰值 .

Setup

我将我的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 += AB 时将 C 矩阵归零 . 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 code

我的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);
}

Assmebly (Intel Syntax)

我的第一次尝试是为这个部分创建一个单独的汇编例程，并从 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"
    );
}

Results

我使用以下标志编译我的代码：

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

Questions

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

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

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

• 为什么我在第一组结果中看到44个Gflops？这高于我计算的12 Gflops峰值，如果我使用单精度计算运行两个核心，这就是我所期望的 .

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