如何通过点积产品达到峰值CPU性能？

matmul 发表于 Dev

Matmul

问题

我一直在研究HPC，特别是使用矩阵乘法作为我的项目（请参阅个人资料中的其他文章）。我在这些方面取得了不错的成绩，但还不够好。我要退后一步，以查看点积计算的效果。

点积与矩阵乘法

点产品比较简单，可以让我测试HPC概念，而无需处理包装和其他相关问题。缓存阻止仍然是一个问题，这构成了我的第二个问题。

算法

将n两个double数组中的对应元素相乘A，B然后求和。一个double组装点积仅仅是一系列的movapd，mulpd，addpd。展开并且设置成在一个聪明的方法，它有可能具有的基团movapd/ mulpd/addpd其上的不同操作xmm的寄存器，并因此独立，优化流水线。当然，事实证明，这无关紧要，因为我的CPU执行不正常。另请注意，重新安排需要剥离最后一次迭代。

其他假设

我不是在编写通用点产品的代码。该代码是针对特定大小的，我不处理附带情况。这只是为了测试HPC概念，并查看我可以获得哪种类型的CPU使用率。

结果

与编译gcc -std=c99 -O2 -m32 -mincoming-stack-boundary=2 -msse3 -mfpmath=sse,387 -masm=intel。我与平常使用的计算机不同。这台计算机的硬件i5 540m可以2.8 GHz * 4 FLOPS/cycle/core = 11.2 GFLOPS/s per core在两步Intel Turbo Boost之后获得（两个内核现在都在运行，因此仅获得2步...如果我关闭一个内核，则可以进行4步Boost）。设置为使用一个线程运行时，32位LINPACK的速度约为9.5 GFLOPS / s。

       N   Total Gflops/s         Residual
     256         5.580521    1.421085e-014
     384         5.734344   -2.842171e-014
     512         5.791168    0.000000e+000
     640         5.821629    0.000000e+000
     768         5.814255    2.842171e-014
     896         5.807132    0.000000e+000
    1024         5.817208   -1.421085e-013
    1152         5.805388    0.000000e+000
    1280         5.830746   -5.684342e-014
    1408         5.881937   -5.684342e-014
    1536         5.872159   -1.705303e-013
    1664         5.881536    5.684342e-014
    1792         5.906261   -2.842171e-013
    1920         5.477966    2.273737e-013
    2048         5.620931    0.000000e+000
    2176         3.998713    1.136868e-013
    2304         3.370095   -3.410605e-013
    2432         3.371386   -3.410605e-013

问题1

我怎么能做得更好呢？我什至无法接近最高表现。我已经将汇编代码优化到了天堂。进一步展开可能会稍微提高一点，但展开次数减少似乎会降低性能。

问题2

当为时n > 2048，您会看到性能下降。这是因为我的L1缓存为32KB，而当n = 2048和A和B的double，他们采取了整个缓存。任何更大的数据都会从内存中流式传输。

我尝试了缓存阻止（未在源代码中显示），但是也许我做错了。谁能提供一些代码或解释如何为缓存阻止点积？

源代码

    #include <stdio.h>
    #include <time.h>
    #include <stdlib.h>
    #include <string.h>
    #include <x86intrin.h>
    #include <math.h>
    #include <omp.h>
    #include <stdint.h>
    #include <windows.h>

    // computes 8 dot products
#define KERNEL(address) \
            "movapd     xmm4, XMMWORD PTR [eax+"#address"]      \n\t" \
            "mulpd      xmm7, XMMWORD PTR [edx+48+"#address"]   \n\t" \
            "addpd      xmm2, xmm6                              \n\t" \
            "movapd     xmm5, XMMWORD PTR [eax+16+"#address"]   \n\t" \
            "mulpd      xmm4, XMMWORD PTR [edx+"#address"]      \n\t" \
            "addpd      xmm3, xmm7                              \n\t" \
            "movapd     xmm6, XMMWORD PTR [eax+96+"#address"]   \n\t" \
            "mulpd      xmm5, XMMWORD PTR [edx+16+"#address"]   \n\t" \
            "addpd      xmm0, xmm4                              \n\t" \
            "movapd     xmm7, XMMWORD PTR [eax+112+"#address"]  \n\t" \
            "mulpd      xmm6, XMMWORD PTR [edx+96+"#address"]   \n\t" \
            "addpd      xmm1, xmm5                              \n\t"

#define PEELED(address) \
            "movapd     xmm4, XMMWORD PTR [eax+"#address"]      \n\t" \
            "mulpd      xmm7, [edx+48+"#address"]               \n\t" \
            "addpd      xmm2, xmm6                  \n\t" \
            "movapd     xmm5, XMMWORD PTR [eax+16+"#address"]   \n\t" \
            "mulpd      xmm4, XMMWORD PTR [edx+"#address"]      \n\t" \
            "addpd      xmm3, xmm7                  \n\t" \
            "mulpd      xmm5, XMMWORD PTR [edx+16+"#address"]   \n\t" \
            "addpd      xmm0, xmm4                  \n\t" \
            "addpd      xmm1, xmm5                  \n\t"

inline double 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(64))) ddot_ref(
    int n,
    const double* restrict A,
    const double* restrict B)
{
    double sum0 = 0.0;
    double sum1 = 0.0;
    double sum2 = 0.0;
    double sum3 = 0.0;
    double sum;
    for(int i = 0; i < n; i+=4) {
        sum0 += *(A + i  ) * *(B + i  );
        sum1 += *(A + i+1) * *(B + i+1);
        sum2 += *(A + i+2) * *(B + i+2);
        sum3 += *(A + i+3) * *(B + i+3);
    }
    sum = sum0+sum1+sum2+sum3;
    return(sum);
}

inline double 
__attribute__ ((gnu_inline))        
__attribute__ ((aligned(64))) ddot_asm
(   int n,
    const double* restrict A,
    const double* restrict B)
{

        double sum;

            __asm__ __volatile__
        (
            "mov        eax, %[A]                   \n\t"
            "mov        edx, %[B]                   \n\t"
            "mov        ecx, %[n]                   \n\t"
            "pxor       xmm0, xmm0                  \n\t"
            "pxor       xmm1, xmm1                  \n\t"
            "pxor       xmm2, xmm2                  \n\t"
            "pxor       xmm3, xmm3                  \n\t"
            "movapd     xmm6, XMMWORD PTR [eax+32]  \n\t"
            "movapd     xmm7, XMMWORD PTR [eax+48]  \n\t"
            "mulpd      xmm6, XMMWORD PTR [edx+32]  \n\t"
            "sar        ecx, 7                      \n\t"
            "sub        ecx, 1                      \n\t" // peel
            "L%=:                                   \n\t"
            KERNEL(64   *   0)
            KERNEL(64   *   1)
            KERNEL(64   *   2)
            KERNEL(64   *   3)
            KERNEL(64   *   4)
            KERNEL(64   *   5)
            KERNEL(64   *   6)
            KERNEL(64   *   7)
            KERNEL(64   *   8)
            KERNEL(64   *   9)
            KERNEL(64   *   10)
            KERNEL(64   *   11)
            KERNEL(64   *   12)
            KERNEL(64   *   13)
            KERNEL(64   *   14)
            KERNEL(64   *   15)
            "lea        eax, [eax+1024]             \n\t"
            "lea        edx, [edx+1024]             \n\t"
            "                                       \n\t"
            "dec        ecx                         \n\t"
            "jnz        L%=                         \n\t" // end loop
            "                                       \n\t"
            KERNEL(64   *   0)
            KERNEL(64   *   1)
            KERNEL(64   *   2)
            KERNEL(64   *   3)
            KERNEL(64   *   4)
            KERNEL(64   *   5)
            KERNEL(64   *   6)
            KERNEL(64   *   7)
            KERNEL(64   *   8)
            KERNEL(64   *   9)
            KERNEL(64   *   10)
            KERNEL(64   *   11)
            KERNEL(64   *   12)
            KERNEL(64   *   13)
            KERNEL(64   *   14)
            PEELED(64   *   15)
            "                                       \n\t"
            "addpd      xmm0, xmm1                  \n\t" // summing result
            "addpd      xmm2, xmm3                  \n\t"
            "addpd      xmm0, xmm2                  \n\t" // cascading add
            "movapd     xmm1, xmm0                  \n\t" // copy xmm0
            "shufpd     xmm1, xmm0, 0x03            \n\t" // shuffle
            "addsd      xmm0, xmm1                  \n\t" // add low qword
            "movsd      %[sum], xmm0                \n\t" // mov low qw to sum
            : // outputs
            [sum]   "=m"    (sum)
            : // inputs
            [A] "m" (A),
            [B] "m" (B), 
            [n] "m" (n)
            : //register clobber
            "memory",
            "eax","ecx","edx","edi",
            "xmm0","xmm1","xmm2","xmm3","xmm4","xmm5","xmm6","xmm7"
            );
        return(sum);
}

int main()
{
    // timers
    LARGE_INTEGER frequency, time1, time2;
    double time3;
    QueryPerformanceFrequency(&frequency);
    // clock_t time1, time2;
    double gflops;

    int nmax = 4096;
    int trials = 1e7;
    double sum, residual;
    FILE *f = fopen("soddot.txt","w+");

    printf("%16s %16s %16s\n","N","Total Gflops/s","Residual");
    fprintf(f,"%16s %16s %16s\n","N","Total Gflops/s","Residual");

    for(int n = 256; n <= nmax; n += 128 ) {
        double* A = NULL;
        double* B = NULL;
        A = _mm_malloc(n*sizeof(*A), 64); if (!A) {printf("A failed\n"); return(1);}
        B = _mm_malloc(n*sizeof(*B), 64); if (!B) {printf("B failed\n"); return(1);}

        srand(time(NULL));

        // create arrays
        for(int i = 0; i < n; ++i) {
            *(A + i) = (double) rand()/RAND_MAX;
            *(B + i) = (double) rand()/RAND_MAX;
        }

        // warmup
        sum = ddot_asm(n,A,B);

        QueryPerformanceCounter(&time1);
        // time1 = clock();
        for (int count = 0; count < trials; count++){
            // sum = ddot_ref(n,A,B);
            sum = ddot_asm(n,A,B);
        }
        QueryPerformanceCounter(&time2);
        time3 = (double)(time2.QuadPart - time1.QuadPart) / frequency.QuadPart;
        // time3 = (double) (clock() - time1)/CLOCKS_PER_SEC;
        gflops = (double) (2.0*n*trials)/time3/1.0e9;
        residual = ddot_ref(n,A,B) - sum;
        printf("%16d %16f %16e\n",n,gflops,residual);
        fprintf(f,"%16d %16f %16e\n",n,gflops,residual);

        _mm_free(A);
        _mm_free(B);
    }
    fclose(f);
    return(0); // successful completion
}

编辑：大会的解释

点积只是两个数字的乘积的重复和：sum += a[i]*b[i]。sum必须0在第一次迭代之前初始化为。向量化后，您一次可以进行2次加和，必须在最后加和：[sum0 sum1] = [a[i] a[i+1]]*[b[i] b[i+1]]，sum = sum0 + sum1。在（英特尔）组装中，这是3个步骤（初始化之后）：

pxor   xmm0, xmm0              // accumulator [sum0 sum1] = [0 0]
movapd xmm1, XMMWORD PTR [eax] // load [a[i] a[i+1]] into xmm1
mulpd  xmm1, XMMWORD PTR [edx] // xmm1 = xmm1 * [b[i] b[i+1]]
addpd  xmm0, xmm1              // xmm0 = xmm0 + xmm1

在这一点上，您没有什么特别的，编译器可以提出。通常，通过将代码展开足够的次数以使用所有xmm可用的寄存器（32位模式下为8个寄存器），可以获得更好的性能。所以，如果你把它打开4倍，让您能够利用所有8个寄存器xmm0通过xmm7。您将有4个累加器和4个寄存器用于存储的结果movapd和addpd。同样，编译器可以提出这个建议。真正思考的部分正在尝试提出一种对代码进行流水线化的方法，即使MOV / MUL / ADD组中的每条指令在不同的寄存器上运行，以便所有3条指令同时执行（通常情况下大多数CPU）。这就是您击败编译器的方式。因此，您必须对4倍展开代码进行图案处理才能做到这一点，这可能需要提前加载向量并剥离第一个或最后一个迭代。这是什么KERNEL(address)是。为了方便起见，我制作了一个4x展开流水线代码的宏。这样，我只需更改即可轻松将其展开为4的倍数address。每个KERNEL计算8个点积。

玻色子

要回答您的总体问题，您无法通过点积获得最佳性能。

问题在于您的CPU可以在每个时钟周期执行一次128位加载，而要进行点积运算，则每个时钟周期需要进行两次128位加载。

但这比大n还要差。您的第二个问题的答案是，点积是受内存限制的，而不是受计算限制的，因此它不能与具有快速核的大n并行化。这在这里更好地解释了为什么对循环进行矢量化没有性能改进。快速内核并行化是一个大问题。我花了一些时间才弄清楚这一点，但学习非常重要。

实际上，几乎没有什么基本算法可以完全受益于快速内核上的并行化。就BLAS算法而言，只有第3级算法（O（n ^ 3））（例如矩阵乘法）才真正受益于并行化。在较慢的内核（例如GPU和Xeon Phi）上，情况会更好，因为内存速度和内核速度之间的差异要小得多。

如果您想找到一种算法，可以使小n接近峰值触发器，请尝试使用例如标量*向量或标量*向量之和。第一种情况应在每个时钟周期进行一次加载，一个多路运算和一个存储，第二种情况应在每个时钟周期进行一次多路，一个加法和一个加载。

我在Knoppix 7.3 32位的Core 2 Duo [email protected]上测试了以下代码。对于标量积，我得到的峰值约为75％，对于标量积之和，我得到的峰值约为75％。标量积的flops /周期为2，标量积之和为4。

编译与 g++ -msse2 -O3 -fopenmp foo.cpp -ffast-math

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <x86intrin.h>

void scalar_product(double * __restrict a, int n) {
    a = (double*)__builtin_assume_aligned (a, 64);
    double k = 3.14159;
    for(int i=0; i<n; i++) {
        a[i] = k*a[i]; 
    }
}

void scalar_product_SSE(double * __restrict a, int n) {
    a = (double*)__builtin_assume_aligned (a, 64);
    __m128d k = _mm_set1_pd(3.14159);    
    for(int i=0; i<n; i+=8) {
        __m128d t1 = _mm_load_pd(&a[i+0]);
        _mm_store_pd(&a[i],_mm_mul_pd(k,t1));
        __m128d t2 = _mm_load_pd(&a[i+2]);
        _mm_store_pd(&a[i+2],_mm_mul_pd(k,t2));
        __m128d t3 = _mm_load_pd(&a[i+4]);
        _mm_store_pd(&a[i+4],_mm_mul_pd(k,t3));
        __m128d t4 = _mm_load_pd(&a[i+6]);
        _mm_store_pd(&a[i+6],_mm_mul_pd(k,t4));
    }
}

double scalar_sum(double * __restrict a, int n) {
    a = (double*)__builtin_assume_aligned (a, 64);
    double sum = 0.0;    
    double k = 3.14159;
    for(int i=0; i<n; i++) {
        sum += k*a[i]; 
    }
    return sum;
}

double scalar_sum_SSE(double * __restrict a, int n) {
    a = (double*)__builtin_assume_aligned (a, 64);
    __m128d sum1 = _mm_setzero_pd();
    __m128d sum2 = _mm_setzero_pd();
    __m128d sum3 = _mm_setzero_pd();
    __m128d sum4 = _mm_setzero_pd();
    __m128d k = _mm_set1_pd(3.14159);   
    for(int i=0; i<n; i+=8) {
        __m128d t1 = _mm_load_pd(&a[i+0]);
        sum1 = _mm_add_pd(_mm_mul_pd(k,t1),sum1);
        __m128d t2 = _mm_load_pd(&a[i+2]);
        sum2 = _mm_add_pd(_mm_mul_pd(k,t2),sum2);
        __m128d t3 = _mm_load_pd(&a[i+4]);
        sum3 = _mm_add_pd(_mm_mul_pd(k,t3),sum3);
        __m128d t4 = _mm_load_pd(&a[i+6]);
        sum4 = _mm_add_pd(_mm_mul_pd(k,t4),sum4);      
    }
    double tmp[8];
    _mm_storeu_pd(&tmp[0],sum1);
    _mm_storeu_pd(&tmp[2],sum2);
    _mm_storeu_pd(&tmp[4],sum3);
    _mm_storeu_pd(&tmp[6],sum4);
    double sum = 0;
    for(int i=0; i<8; i++) sum+=tmp[i];
    return sum;
}

int main() {
    //_MM_SET_FLUSH_ZERO_MODE(_MM_FLUSH_ZERO_ON);
    //_mm_setcsr(_mm_getcsr() | 0x8040);
    double dtime, peak, flops, sum;
    int repeat = 1<<18;
    const int n = 2048;
    double *a = (double*)_mm_malloc(sizeof(double)*n,64);
    double *b = (double*)_mm_malloc(sizeof(double)*n,64);
    for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;

    dtime = omp_get_wtime();
    for(int r=0; r<repeat; r++) {
        scalar_product_SSE(a,n);
    }
    dtime = omp_get_wtime() - dtime;
    peak = 2*2.67;
    flops = 1.0*n/dtime*1E-9*repeat;
    printf("time %f, %f, %f\n", dtime,flops, flops/peak);

    //for(int i=0; i<n; i++) a[i] = 1.0*rand()/RAND_MAX;
    //sum = 0.0;    
    dtime = omp_get_wtime();
    for(int r=0; r<repeat; r++) {
        scalar_sum_SSE(a,n);
    }
    dtime = omp_get_wtime() - dtime;
    peak = 2*2*2.67;
    flops = 2.0*n/dtime*1E-9*repeat;
    printf("time %f, %f, %f\n", dtime,flops, flops/peak);
    //printf("sum %f\n", sum);

}

本文收集自互联网，转载请注明来源。

如有侵权，请联系[email protected] 删除。

编辑于2020-11-26

我来说两句

0条评论

登录后参与评论

上一篇：Swift：将CInt转换为Int？

来自分类Dev

Related 相关文章

文章

如何通过点积产品达到峰值CPU性能？

如何通过点积产品达到峰值CPU性能？

当CPU通过IIS高级CPU设置达到峰值时杀死W3WP

如何减少PowerShell CPU峰值

如何诊断 Ubuntu CPU 峰值/IO 等待？

核心矩阵Clojure点产品性能

使用self.delegate子类化UITextField导致应用程序冻结，CPU峰值达到100％

如何求和列的成对点积

硬盘性能峰值

CPU温度是否可以导致CPU在任务管理器中显示为峰值达到100％？

Akka如何达到如此高的性能？

从DRAM和L1缓存获取CPU峰值FP性能和延迟以及带宽

RaspberryPi cpu温度峰值

当游戏中的CPU使用率达到峰值时，为什么我的NAudio音乐会冻结？

带有SSE说明的点产品性能

Linux +如何验证CPU性能

多个向量对的块和点积：如何完成？

如何矢量化轨迹段之间的点积

如何找到两行的点积（Opencv）

如何使用numpy的einsum取子数组的点积？

如何在numpy的矩阵向量中制作点积？

TensorFlow：如何对向量和张量进行点积？

如何在Caffe的矩阵之间做点积？

如何使用点积计算权重的梯度下降成本？

Numpy：如何在张量上正确执行点积

点积稀疏矩阵

Matlab：集成点积

fft的点积

向量的点积

numpy点积逐步

熊猫矢量的点积