如何获得基本向量/矩阵运算的理论/Linpack失败性能？

Question

问题很简单。如何进一步优化我的代码，因为基本的矩阵操作对于我的计算来说是非常重要和常见的。BLAS和LAPACK运算在线性代数中是很好的，但它们都不提供基本的元素加法/乘运算(Hadamard)。理论性能可能比较困难，但Linpack性能或60%~80%的Linpack性能应该是可以实现的。(我只能做12%，如果我使用乘加，那么只有25%)

参考资料

理论性能:8259 u有4个核* 3.8GHz * 16触发器= 240 GFlops

Linpack性能:8259 u可运行140~160 GFlops双精度操作。

平台: Macbook Pro 2018，蒙特里

CPU: i5-8259 u，4c8t

RAM: 8GB

CC: gcc 11.3.0

-mavx2 -mfma -fopenmp -O3

，这是我的尝试，

触发器的计算如下：

double time = stop - start;
double ops = 1.0 * Nx * Ny * iterNum; //2.0 for complex numbers
double flops = ops / time;
double gFlops = flops / 1E9;

下面是我运行代码时的一些结果。真实的和复杂的结果几乎是一样的。只显示实际结果(大致)：

//Nx = Ny = 2048, iterNum = 10000
//Typical matrix size and iteration depth for my calculation
threads =    1:  1 GFlops
threads =    2:  2 GFlops
threads =    4:  3 GFlops
threads =    8:  4 GFlops
threads =   16:  9 GFlops
threads =   32: 11 GFlops
threads =   64: 15 GFlops
threads =  128: 18 GFlops
threads =  256: 19 GFlops
threads =  512: 21 GFlops
threads = 1024: 20 GFlops
threads = 2048: 40 GFlops // wrong answer

为了便于堆上的大型矩阵和与mathGL的集成，将矩阵作为由Nx * Ny元素组成的向量按行级联。

// for real numbers
  x = (double *)_mm_malloc(Nx * Ny * sizeof(double), 32);
  y = (double *)_mm_malloc(Nx * Ny * sizeof(double), 32);
  z = (double *)_mm_malloc(Nx * Ny * sizeof(double), 32);
sum = (double *)_mm_malloc(Nx * Ny * sizeof(double), 32);
// for complex numbers
  x = (double *)_mm_malloc(Nx * Ny * sizeof(double complex), 32);
  y = (double *)_mm_malloc(Nx * Ny * sizeof(double complex), 32);
  z = (double *)_mm_malloc(Nx * Ny * sizeof(double complex), 32);
sum = (double *)_mm_malloc(Nx * Ny * sizeof(double complex), 32);

同时用openmp进行加成。

    double start = omp_get_wtime();

    #pragma omp parallel private(shift)
    {
        for (int tds = omp_get_thread_num(); tds < threads; tds = tds + threads)
        {
            shift = Nx * Ny / threads * tds;
            for (int i = 0; i < iterNum; i++)
            {
                AddComplex(sum+shift, sum+shift, z+shift, Nx/threads, Ny);
            }
        }
    }

    double stop = omp_get_wtime();

我使用AVX本质"immintrin.h“编写了显式向量化代码。

//real matrix addition
void AddReal(double *summation, const double *summand, const double *addend, int Nx, int Ny)
{
    int nBlock = Nx * Ny / realPackSize;
    int nRem = Nx * Ny % realPackSize;

    register __m256d packSummand, packAddend, packSum;
    const double *px = summand;
    const double *py = addend;
    double *pSum = summation;
    
    for (int i = 0; i < nBlock; i++)
    {
        packSummand = _mm256_load_pd(px);
        packAddend = _mm256_load_pd(py);
        packSum = _mm256_add_pd(packSummand, packAddend);
        _mm256_store_pd(pSum, packSum);
        px = px + realPackSize;
        py = py + realPackSize;
        pSum = pSum + realPackSize;
    }
    for (int i = 0; i < nRem; i++)
    {
        pSum[i] = px[i] + py[i];
    }

    px = NULL;
    py = NULL;
    pSum = NULL;

    return;
}
//Complex matrix addition
void AddComplex(double complex *summation, const double complex *summand, const double complex *addend, int Nx, int Ny)
{
    int nBlock = Nx * Ny / complexPackSize;
    int nRem = Nx * Ny % complexPackSize;

    register __m256d packSummand, packAddend, packSum;
    const double complex *px = summand;
    const double complex *py = addend;
    double complex *pSum = summation;

    for (int i = 0; i < nBlock; i++)
    {
        packSummand = _mm256_load_pd(px);
        packAddend = _mm256_load_pd(py);
        packSum = _mm256_add_pd(packSummand, packAddend);
        _mm256_store_pd(pSum, packSum);
        px = px + complexPackSize;
        py = py + complexPackSize;
        pSum = pSum + complexPackSize;
    }
    for (int i = 0; i < nRem; i++)
    {
        pSum[i] = px[i] + py[i];
    }
    
    px = NULL;
    py = NULL;
    pSum = NULL;

    return;
}

Answer 1

1级(如点积)和等级2(例如。已知BLAS函数不缩放(特别是1级BLAS函数)，而不是3级(例如。矩阵-乘法)。实际上，它们通常是内存绑定的：读取/写入的数据量是O(n)，而浮点操作的数量也是O(n)。对于级别3的BLAS来说，情况并非如此，后者通常是明确的计算范围。

理论性能可能比较困难，但Linpack性能或60%~80%的Linpack性能应该是可以实现的。

如果计算是内存绑定，那么，不，这是不可能的。在几乎所有的机器上，Linpack通常都是明确的计算界。在过去的几十年里，内存的速度并没有处理器的速度增长的那么快。这被称为内存墙(几十年前制定的，现在仍然如此)。

，这是我运行代码时的一些结果。

使用1024个线程，而不是在一个有4个核心和8个线程的移动处理器上使用512，计算速度更快，这让我认为在某个地方存在着一个巨大的问题。最大值应该达到8个线程，否则，这意味着计算显然是低效的。事实上，运行比硬件线程更多的线程会导致操作系统调度程序进行昂贵的上下文切换(更高的开销)。最后，您的处理器一次运行超过8个任务。有两种可能性：

，

• ，时间不对(所提供的代码对我来说很好)
• ，程序是伪造的
• ，计算显示超线性速度(可能是由于缓存)

。

我使用AVX "immintrin.h“编写了显式向量化代码。

热循环包含两个加载、一个存储、一个添加和少量增加整数的指令。您的处理器可以完成2次加载和1次存储，因此在假设nBlock足够大的情况下，SIMD部分可以在一个吞吐量周期内完成(尽管延迟可能要大得多)。

您的处理器每个周期可以进行2次添加，这样一半的吞吐量就会丢失。但是，如果加载/写入是强制性的，则不能编写比这更快的东西。

如果complexPackSize小于SIMD车道，那么我认为由于与过去的迭代重叠，处理器必须进行复杂的操作，这肯定会大大降低循环的运行效率(带循环的依赖项将使循环延迟界变得非常低效)。如果complexPackSize比缓存行大得多，那么预取可能是一个问题。

您的处理器不能同时执行太多指令。增量指令和循环检查导致5条指令被执行，这消耗了至少一个周期。这再次将吞吐量降低了2倍，因此不超过理论性能的25%。通过展开循环，可以稍微改进这一点。展开也可能改善执行，因为_mm256_add_pd指令有相当高的延迟。应该记住，SIMD指令对于吞吐量来说是很好的，但对于延迟则不是。因此，当延迟不是问题时，SIMD代码应该是快速的。

请注意，在使用_mm256_store_pd时，写分配缓存策略会导致读取数据，增加从RAM传输的数据量，并降低观察到的吞吐量。可以使用_mm256_stream_pd来避免这种影响，但只有在数据没有在缓存后读取或数据不适合缓存时才能快速读取。它还要求数据必须是对齐的。事实上，_mm256_store_pd也要求这样做，如果不是这样的话，它肯定会导致一个沉默的错误。同样适用于_mm256_load_pd：_mm256_loadu_pd应该用于未对齐的数据。我不确定读取的数据总是对齐的。如果complexPackSize是一个可以被32除的二次方以及shift，那应该是很好的。但是，我对shift的情况非常怀疑，特别是在大量线程中。我也发现非常可疑地使用一个恒定的complexPackSize，而SIMD车道有一个固定的大小。你检查过所有病例的结果了吗？