SIMD

SIMD

关注
发信
关注(28)粉丝(399)

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

问题
我正在研究高性能矩阵乘法算法，如openblas或gotoblas，我正试图重现一些结果。这个问题涉及矩阵乘法算法的内核。具体地说，我正在研究计算，其中C += AB和A是类型为B的2x2矩阵，在CPU的峰值速度下。有两种方法可以做到这一点。一种方法是使用simd指令。第二种方法是使用simd寄存器直接在程序集中编码。
到目前为止我所看到的
所有相关的论文，课程网页，许多关于这个主题的问答（太多了，无法列出），我已经在我的电脑上编译了openblas，浏览了openblas，gotoblas，blis源代码，agner的手册。
硬件
我的CPU是Intel i5-540M。您可以在cpu-world.com上找到相关的CPuid信息。微体系结构是nehalem（westmere），因此理论上它可以计算每个核心周期4个双精度触发器。我将只使用一个内核（没有openmp），因此随着超线程关闭和4步英特尔涡轮增压，我应该看到一个峰值double。作为参考，在两个内核都在峰值运行的情况下，“英特尔涡轮增压”提供了2步加速，理论峰值应该是( 2.533 Ghz + 4*0.133 Ghz ) * ( 4 DP flops/core/cycle ) * ( 1 core ) = 12.27 DP Gflops。
安装程序
我将2x2矩阵声明为22.4 DP Gflops并使用随机条目初始化它们，如下面的代码片段所示。

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的估计值，我调用每个乘法函数一次来预热，然后在double循环中执行for次，确保每次计算maxiter时C矩阵都归零。C += AB循环放在两个for语句中，用于估计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位向量，因此我可以在每个向量中拟合2clock()s。这是我在内核中进行2x2矩阵乘法的主要原因。simd代码一次计算一整行double。
    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（英特尔语法）
我的第一个尝试是为这个部分创建一个单独的程序集例程，并从C例程调用它。但是，它非常慢，因为我不能内联main函数。我将程序集编写为内联程序集，如下所示。它与extern产生的相同。根据我对Nehalem微体系结构图的理解，该处理器可以并行执行gcc -S -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel、SSE ADD和SSE MUL，这解释了SSE MOV、MUL、ADD指令的交错。您会注意到上面的simd指令顺序不同，因为我对agner fog的手册有不同的理解。不过，MOV很聪明，上面的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"
    );
}

结果
我使用以下标志编译代码：
gcc -std=c99 -O3 -msse3 -ffast-math -march=nocona -mtune=nocona -funroll-all-loops -fomit-frame-pointer -masm=intel

gcc的结果如下：
********** 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版本不与maxiter = 1000000000内联，结果是：
********** 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与“noinline”simd性能相同的结果中可以明显看出。我能找到的唯一解释是《阿涅尔雾》第2卷第6页：
编译的代码可能比汇编代码快，因为编译器可以
程序间优化和整个程序优化。大会
程序员通常必须使用定义良好的调用来创建定义良好的函数
遵循所有调用约定以使代码可测试的接口
可核实。这就阻止了编译器使用的许多优化方法，例如
作为函数内联、寄存器分配、常量传播、公共子表达式
跨功能消除、跨功能调度等
通过使用具有内在功能的C++代码，可以获得优点。
装配代码。
但两个版本的汇编程序输出完全相同。
为什么我在第一组结果中看到44个gflop？这比我计算的12 gflops峰值要高得多，如果我用单精度计算运行两个内核，这就是我所期望的。
编辑1
评论说可能存在死代码消除，我可以确认，这是发生在simd指令。__attribute__ ((noinline))输出显示simd的-S循环仅为0for矩阵。我可以用C关闭编译器优化来禁用它。在这种情况下，simd的运行速度是asm的3倍，但asm仍然以完全相同的速度运行。现在的范数也不是零，但是在10^-16还是可以的。我还看到内联asm版本是用-O0和APP标记内联的，但它也在NO_APP循环中展开了8次。我认为多次展开将严重影响性能，因为我通常展开循环4次。以我的经验来看，任何事情都会降低性能。

最佳答案

gcc正在使用intrinsic优化内联函数，mult2by2B，因为

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

如果没有那条线，科利鲁在电脑上花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);
}

关于c - 优化的2x2矩阵乘法:慢速组装与快速SIMD，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/23790842/

10-11 18:36
Powered by LMLPHP ©2024 SIMD 0.032273