c++ - 如何使用 Vector Class Library 进行 AVX 矢量化以及 openmp #pragma omp parallel 进行缩减？

Question

我正在使用 OpenMP 并行化循环，即在内部使用 AVX-512 和 Agner Fog's VCL Vector Class Library.

这是代码:

double HarmonicSeries(const unsigned long long int N) {
  unsigned long long int i;
  Vec8d divV(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);
  Vec8d sumV(0.0);
  const Vec8d addV(8.0);
  const Vec8d oneV(1.0);

  #pragma omp parallel for reduction(+:sumV,divV)
  for(i=0; i<N; ++i) {
    sumV += oneV / divV;
    divV += addV;
  }
  return horizontal_add(sumV);
}

当尝试编译上面的代码时，我得到了

g++ -Wall -Wextra -O3 -g -I include -fopenmp -m64 -mavx2 -mfma -std=c++17 -o harmonic_series harmonic_series.cpp
harmonic_series.cpp:87:40: error: user defined reduction not found for ‘sumV’
   87 |   #pragma omp parallel for reduction(+:sumV,divV)
      |                                        ^~~~
harmonic_series.cpp:87:45: error: user defined reduction not found for ‘divV’
   87 |   #pragma omp parallel for reduction(+:sumV,divV)

有关如何解决此问题并为 Vec8d 类提供用户定义的缩减有任何提示吗？它只是由 VCL 类定义的加号运算符，但我找不到任何如何编写此代码的示例。

非常感谢您的帮助!

Answer 1

这是最终的解决方案。它使用自动归约，并通过仅在每个循环的第一次迭代中计算 divV = startdivV + i * addV ，然后使用 divV += addV for 来避免 u64->fp 转换所有其他迭代。计算前 9.6e10 元素之和的运行时为 {real 9s, user 1m46s}，在 Intel Core i7-10850H CPU 上具有 12 个线程 - 与 manual reduction. 相同。

double HarmonicSeries(const unsigned long long int N) {
  unsigned long long int i;
  Vec8d divV(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0);
  Vec8d sumV(0.0);
  const Vec8d addV(8.0);
  const Vec8d oneV(1.0);
  const Vec8d startdivV = divV;
  bool first_loop = true;
  #pragma omp declare reduction( + : Vec8d : omp_out = omp_out + omp_in ) initializer (omp_priv=omp_orig)
//It's important to mark "first_loop" variable as firstprivate so that each private copy gets initialized.
  #pragma omp parallel for firstprivate(first_loop) lastprivate(divV) reduction(+:sumV)
  for(i=0; i<N; ++i) {
    if (first_loop) {
      divV = startdivV + i * addV;
      first_loop = false;
    } else {
      divV += addV;
    }
    sumV += oneV / divV;
  }
  return horizontal_add(sumV);
}