c++ - 在 openmp 中创建一个 for 循环，与矩阵/vector 操作并行

Question

我有这个代码:

scalar State::add(const int N, const int M,                                                                                                                                                                                                                                               
                        vector<scalar>& flmn,                                                                                                                                                                                                                                             
                        vector<scalar>& BSum,                                                                                                                                                                                                                                             
                        const vector<scalar>& prev_flm,                                                                                                                                                                                                                                   
                        const vector<scalar>& prev_bigsum,                                                                                                                                                                                                                                
                        const vector<scalar>& Qratio,                                                                                                                                                                                                                                     
                        const int test)                                                                                                                                                                                                                                                   
{                                                                                                                                                                                                                                                                                         
  scalar c=1;                                                                                                                                                                                                                                                                          
  #pragma omp parallel for                                                                                                                                                                                                                                                                
  for(int i=1;i<=M;i++)                                                                                                                                                                                                                                                                   
   {                                                                                                                                                                                                                                                                                      
     flmn.at(i-1) = Qratio.at(i-1)*k1+k2;                                                                                                                                                                                                                                                 
     BSum.at(i-1) = someconstant + somepublicvector.at(1)*flmn.at(i-1);                                                                                                                                                                                                                   
     c *= BSum.at(i-1);                                                                                                                                                                                                                                                                   
   }                                                                                                                                                                                                                                                                                      
   return c;                                                                                                                                                                                                                                                                              
}                                                                                                                                                                                                                                                                                         

最后我返回变量 c。当使用这个:“#pragma omp parallel for”时，它肯定不会给我一致的答案，因为迭代之间总是存在重叠。我想知道如何在 openmp 中并行化这种矩阵或 vector 操作的组合，而且我会得到一致的结果，因为这里显然存在竞争条件问题？

Answer 1

for (int i = 1; i <= M; i++) {
    flmn.at(i - 1) = Qratio.at(i - 1) * k1 + k2;
    BSum.at(i - 1) = someconstant + somepublicvector.at(1) * flmn.at(i - 1);
    c *= BSum.at(i - 1);
}

一些注意事项:

1. 除非您确实需要异常安全索引，否则不要使用 std::vector::at。

2. 您对每个 vector 使用相同的索引，因此从 i=0 而不是 Fortran 风格的 i=1 开始。

3. M 是否与所用 vector 的大小不同(即它是一个子集)？如果不是，则无需指定。

---

一个可能的 OpenMP 实现可以是

scalar c{1.0};

#pragma omp parallel
{
    const std::size_t nthreads = omp_get_num_threads();
    const std::size_t chunk_size = M / nthreads;    // WARNING: non-even division case left to user
    const std::size_t tid = omp_get_thread_num();

    #pragma omp for reduction(*:c)
    for (std::size_t j = 0; j < chunk_size; j++) {
        const std::size_t i = j + tid * chunk_size;
        flmn[i] = Qratio[i] * k1 + k2;
        BSum[i] = someconstant + somepublicvector[1] * flmn[i];
        c *= BSum[i];
    }
}

请注意，我假设 nthreads 平均划分 M。如果不是，则需要单独处理此案例。如果您使用的是 OpenMP 4.0，那么我建议使用 simd 指令，因为前两行都是 saxpy 操作并且可以从矢量化中获益。为获得最佳性能，请确保 chunk_size 是 CPU 缓存行大小的倍数。