OpenMP/C++：并行 for 循环，之后进行归约 - 最佳实践？

Question

给定以下代码...

for (size_t i = 0; i < clusters.size(); ++i)
{
    const std::set<int>& cluster = clusters[i];
    // ... expensive calculations ...
    for (int j : cluster)
        velocity[j] += f(j); 
}

...我想在多个 CPUs/cores 上运行。函数 f 不使用 velocity.

在第一个 for 循环之前的一个简单的 #pragma omp parallel for 将产生 unpredictable/wrong 结果，因为 std::vector<T> velocity 在内部循环中被修改。多个线程可以同时访问和（尝试）修改 velocity 的相同元素。

我认为第一个解决方案是在velocity[j] += f(j);操作之前写#pragma omp atomic。这给了我一个编译错误（可能与 Eigen::Vector3d 或 velocity 类型的元素是 class 成员有关）。另外，我读到原子操作非常慢与每个线程都有一个私有变量并在最后进行减少相比。所以这就是我想做的，我想。

我想到了这个：

#pragma omp parallel
{
    // these variables are local to each thread
    std::vector<Eigen::Vector3d> velocity_local(velocity.size());
    std::fill(velocity_local.begin(), velocity_local.end(), Eigen::Vector3d(0,0,0));

    #pragma omp for
    for (size_t i = 0; i < clusters.size(); ++i)
    {
        const std::set<int>& cluster = clusters[i];
        // ... expensive calculations ...
        for (int j : cluster)
            velocity_local[j] += f(j); // save results from the previous calculations
    } 

    // now each thread can save its results to the global variable
    #pragma omp critical
    {
        for (size_t i = 0; i < velocity_local.size(); ++i)
            velocity[i] += velocity_local[i];
    }
}

这是一个好的解决方案吗？它是最佳解决方案吗？（甚至正确？）

进一步思考：使用 reduce 子句（而不是 critical 部分）会引发编译器错误。我认为这是因为 velocity 是 class 成员。

找了一个问题类似的问题，this问题貌似差不多。但我认为我的情况可能不同，因为最后一步包含一个 for 循环。这是否是最佳方法的问题仍然存在。

编辑： 根据评论要求：reduction 子句...

    #pragma omp parallel reduction(+:velocity)
    for (omp_int i = 0; i < velocity_local.size(); ++i)
        velocity[i] += velocity_local[i];

...抛出以下错误：

错误 C3028：'ShapeMatching::velocity'：数据共享子句中只能使用变量或静态数据成员

（与 g++ 类似的错误）

Answer 1

您正在进行数组缩减。我已经多次描述过这一点（例如 reducing an array in openmp and fill histograms array reduction in parallel with openmp without using a critical section）。您可以在有和没有关键部分的情况下执行此操作。

您已经在关键部分（在您最近的编辑中）正确地做到了这一点，所以让我描述一下如何在没有关键部分的情况下做到这一点。

std::vector<Eigen::Vector3d> velocitya;
#pragma omp parallel
{
    const int nthreads = omp_get_num_threads();
    const int ithread = omp_get_thread_num();
    const int vsize = velocity.size();

    #pragma omp single
    velocitya.resize(vsize*nthreads);
    std::fill(velocitya.begin()+vsize*ithread, velocitya.begin()+vsize*(ithread+1), 
              Eigen::Vector3d(0,0,0));

    #pragma omp for schedule(static)
    for (size_t i = 0; i < clusters.size(); i++) {
        const std::set<int>& cluster = clusters[i];
        // ... expensive calculations ...
        for (int j : cluster) velocitya[ithread*vsize+j] += f(j);
    } 

    #pragma omp for schedule(static)
    for(int i=0; i<vsize; i++) {
        for(int t=0; t<nthreads; t++) {
            velocity[i] += velocitya[vsize*t + i];
        }
    }
}

此方法需要额外的 care/tuning，因为我没有进行过虚假共享。

至于哪种方法更好，你自己去试验吧。

OpenMP/C++：并行 for 循环，之后进行归约 - 最佳实践？

OpenMP/C++: Parallel for loop with reduction afterwards - best practice?

c++

parallel-processing

multithreading

openmp