在矩阵向量乘法中使用 OpenMP“for simd”？

如何解决在矩阵向量乘法中使用 OpenMP“for simd”？

我目前正在尝试通过将 #pragma omp for 与 #pragma omp simd 结合来使我的矩阵向量乘法函数与 BLAS 相媲美，但与我仅使用为构造。如何使用 OpenMP 的 SIMD 构造正确矢量化内循环？

vector dot(const matrix& A,const vector& x)
{
  assert(A.shape(1) == x.size());

  vector y = xt::zeros<double>({A.shape(0)});

  int i,j;
#pragma omp parallel shared(A,x,y) private(i,j)
  {
#pragma omp for // schedule(static)
    for (i = 0; i < y.size(); i++) { // row major
#pragma omp simd
      for (j = 0; j < x.size(); j++) {
        y(i) += A(i,j) * x(j);
      }
    }
  }

  return y;
}

解决方法

您的指令不正确，因为会引入竞争条件（在 y(i) 上）。在这种情况下，您应该使用归约。下面是一个例子：

vector dot(const matrix& A,const vector& x)
{
  assert(A.shape(1) == x.size());

  vector y = xt::zeros<double>({A.shape(0)});

  int i,j;

  #pragma omp parallel shared(A,x,y) private(i,j)
  {
    #pragma omp for // schedule(static)
    for (i = 0; i < y.size(); i++) { // row major
      decltype(y(0)) sum = 0;

      #pragma omp simd reduction(+:sum)
      for (j = 0; j < x.size(); j++) {
        sum += A(i,j) * x(j);
      }

      y(i) += sum;
    }
  }

  return y;
}

请注意，可能不需要更快，因为某些编译器能够自动矢量化代码（例如 ICC）。 GCC 和 Clang 经常无法自动执行（高级）SIMD 缩减，这样的指令对他们有点帮助。您可以检查汇编代码以检查代码如何矢量化或启用矢量化报告（请参阅 here 以了解 GCC）。