时间戳

首页 > 解决方案 > 在矩阵向量乘法中使用 OpenMP“for simd”？

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

问题描述

我目前正在尝试通过与 组合使我的矩阵向量乘法函数与 BLAS 进行比较#pragma omp for，#pragma omp simd但与仅使用 for 构造相比，它没有得到任何加速改进。如何使用 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;
}

标签： parallel-processingopenmpvectorizationsimdxtensor

解决方案

您的指令不正确，因为会引入竞争条件（on 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 缩减，这样的指令可以帮助他们一点。您可以检查汇编代码以检查代码是如何矢量化的或启用矢量化报告（有关 GCC，请参见此处）。

推荐阅读