如何以完美的答案同时使用 openmp 和 AVX2？

Question

我使用 OpenMP 和 AVX2 编写了 Matrix-Vector 产品程序。

但是，由于 OpenMP，我得到了错误的答案。 真正的答案是数组 c 的所有值都会变成 100。

我的答案是 98、99 和 100 的混合。

实际代码如下。

我用 -fopenmp、-mavx、-mfma 编译了 Clang。

#include "stdio.h"
#include "math.h"
#include "stdlib.h"
#include "omp.h"
#include "x86intrin.h"

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    int k;
#pragma omp parallel
    {
        __m256d va,vb,vc;
        int i;
#pragma omp for private(i, va, vb, vc) schedule(static)
        for (k = 0; k < l; k++) {
            vb = _mm256_broadcast_sd(&b[k]);
            for (i = 0; i < m; i+=4) {
                va = _mm256_loadu_pd(&a[m*k+i]);
                vc = _mm256_loadu_pd(&c[i]);

                vc = _mm256_fmadd_pd(vc, va, vb);

                _mm256_storeu_pd( &c[i], vc );
            }
        }
    }
}
int main(int argc, char* argv[]) {

    // set variables
    int m;
    double* a;
    double* b;
    double* c;
    int i;

    m=100;
    // main program

    // set vector or matrix
    a=(double *)malloc(sizeof(double) * m*m);
    b=(double *)malloc(sizeof(double) * m*1);
    c=(double *)malloc(sizeof(double) * m*1);
    //preset
    for (i=0;i<m;i++) {
        a[i]=1;
        b[i]=1;
        c[i]=0.0;
    }
    for (i=m;i<m*m;i++) {
        a[i]=1;
    }

    mv(a, b, c, m, 1, m);

    for (i=0;i<m;i++) {
        printf("%e\n", c[i]);
    }
    free(a);
    free(b);
    free(c);
    return 0;
}

我知道关键部分会有所帮助。但是临界区很慢。

那么，我该如何解决这个问题呢？

Answer 1

你想要的基本操作是

c[i] = a[i,k]*b[k]

如果你使用row-major order storage，这就变成了

c[i] = a[i*l + k]*b[k]

如果您使用以列为主的顺序存储，这将变为

c[i] = a[k*m + i]*b[k]

对于行主要顺序，您可以像这样并行化

#pragma omp parallel for
for(int i=0; i<m; i++) {
  for(int k=0; k<l; k++) {
    c[i] += a[i*l+k]*b[k];
  }
}

对于列主顺序，您可以像这样并行化

#pragma omp parallel
for(int k=0; k<l; k++) {
  #pragma omp for
  for(int i=0; i<m; i++) {
    c[i] += a[k*m+i]*b[k];
  }
}

矩阵向量操作是二级操作，是内存带宽限制操作。 1 级和 2 级操作不随内核数量等因素扩展。只有 3 级操作（例如密集矩阵乘法）可以扩展 https://en.wikipedia.org/wiki/Basic_Linear_Algebra_Subprograms#Level_3。

Answer 2

问题不在于您的 AVX 内部函数，让我们看一下没有内部函数的代码：

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    #pragma omp parallel for schedule(static)
    for (int k = 0; k < l; k++) {
        double xb = b[k];
        for (int i = 0; i < m; i++) {
            double xa = a[m*k+i];
            double xc = c[i];
            xc = xc + xa * xb;
            c[i] = xc;
        }
    }
}

注意：您的私有声明在技术上是正确且多余的，因为在并行循环内声明了，但如果您尽可能在本地声明变量，则更容易推理代码。

您的代码的竞争条件是c[i] - 多个线程尝试更新。现在，即使您可以通过原子更新来保护它，性能也会很糟糕：不仅因为保护，还因为 c[i] 的数据必须在不同内核的缓存之间不断移动。

您可以做的一件事是在c 上使用数组缩减。这会为每个线程创建一个 c 的私有副本，并在最后合并：

void mv(double *a,double *b,double *c, int m, int n, int l)
{
    #pragma omp parallel for schedule(static) reduction(+:c[:m])
    for (int k = 0; k < l; k++) {
        for (int i = 0; i < m; i++) {
            c[i] += a[m*k+i] * b[k];
        }
    }
}

只要两个m-vectors 适合你的缓存，这应该是相当有效的，但由于线程管理开销，你仍然可能会得到很多开销。最终，您将受到内存带宽的限制，因为在向量矩阵乘法中，从 a 读取的每个元素只有一次计算。

无论如何，您当然可以交换i 和k 循环并保存减少，但是您在a 上的内存访问模式将是低效的（跨步） - 所以你应该block 循环以避免那个。

现在，如果您查看output of any modern compiler，它将自行生成 SIMD 代码。当然，如果您愿意，您可以应用自己的 SIMD 内在函数。但是，如果 m 不能被 4 整除（原始版本中没有），请确保正确处理边缘情况。

归根结底，如果您真的想要性能，请使用 BLAS 库（例如 MKL）中的函数。如果您想尝试优化，有很多机会深​​入细节。