OpenMP 减少容器元素

Question

我有一个嵌套循环，外部迭代很少，内部迭代很多。在内部循环中，我需要计算一个总和，所以我想使用 OpenMP 缩减。外部循环位于容器上，因此减少应该发生在该容器的元素上。 这是一个最小的人为示例：

#include <omp.h>
#include <vector>
#include <iostream>

int main(){
    constexpr int n { 128 };

    std::vector<int> vec (4, 0);
    for (unsigned int i {0}; i<vec.size(); ++i){

        /* this does not work */
        //#pragma omp parallel for reduction (+:vec[i])
        //for (int j=0; j<n; ++j)
        //  vec[i] +=j;

        /* this works */
        int* val { &vec[0] };
        #pragma omp parallel for reduction (+:val[i])
        for (int j=0; j<n; ++j)
            val[i] +=j;

        /* this is allowed, but looks very wrong. Produces wrong results
         * for std::vector, but on an Eigen type, it worked. */
        #pragma omp parallel for reduction (+:val[i])
        for (int j=0; j<n; ++j)
            vec[i] +=j;
    }
    for (unsigned int i=0; i<vec.size(); ++i) std::cout << vec[i] << " ";
    std::cout << "\n";

    return 0;
}

问题是，如果我将缩减子句写为(+:vec[i])，我会收到错误‘vec’ does not have pointer or array type，它的描述性足以找到解决方法。但是，这意味着我必须引入一个新变量并在一定程度上更改代码逻辑，而且我发现看代码应该做什么不太明显。

我的主要问题是，是否有更好/更清洁/更标准的方法来编写容器元素的缩减。

我也想知道上面代码中显示的第三种方式有点为什么以及如何工作。我实际上正在使用Eigen 库，该变体在其容器上似乎工作得很好（虽然还没有进行广泛的测试），但在std::vector 上，它产生的结果介于零和实际结果之间（8128 )。我认为它应该可以工作，因为vec[i] 和val[i] 应该 都评估为取消引用相同的地址。但很可惜，显然不是。

我正在使用 OpenMP 4.5 和 gcc 9.3.0。

Answer 1

我将分三个部分回答你的问题：

1.在上面的示例中，使用 std::vec 执行 OpenMP 缩减的最佳方法是什么？

i) 使用您的方法，即创建一个指针int* val { &vec[0] };

ii) 声明一个新的共享变量，如 @1201ProgramAlarm 已回答。

iii) 声明 user defined reduction （这在您的简单情况下并不真正适用，但请参阅下面的 3. 以获得更有效的模式）。

2.为什么第三个循环不起作用，为什么它与 Eigen 一起起作用？

就像前面的答案所说，您告诉 OpenMP 对内存地址 X 执行归约和，但您正在对内存地址 Y 执行加法，这意味着归约声明被忽略并且您的加法受制于通常的线程竞争条件。

您并没有真正为您的 Eigen 冒险提供太多细节，但这里有一些可能的解释：

i) 你并没有真正使用多线程（检查n = Eigen::nbThreads( )）

ii) 您没有禁用 Eigen 自己的并行性，这可能会破坏您自己对 OpenMP 的使用，例如EIGEN_DONT_PARALLELIZE 编译器指令。

iii) 存在竞争条件，但您没有看到它，因为 Eigen 操作需要更长的时间，您使用的线程数量较少，并且只写入少量值 => 较少出现的线程干扰每个其他会产生错误的结果。

3.我应该如何使用 OpenMP 并行化这个场景（技术上不是你明确提出的问题）？

您应该同时并行化两者，而不是只并行化内部循环。您拥有的序列号越少越好。在这种情况下，每个线程都有自己的 vec 向量的私有副本，在所有元素都被各自的线程求和后，该向量会减少。此解决方案最适合您提供的示例，但如果您使用非常大的向量和非常多的线程（或 RAM 非常有限），则可能会遇到 RAM 问题。

#pragma omp parallel for collapse(2) reduction(vsum : vec)
for (unsigned int i {0}; i<vec.size(); ++i){
    for (int j = 0; j < n; ++j) {
        vec[i] += j;
    }
}

其中 vsum 是用户定义的归约，即

#pragma omp declare reduction(vsum : std::vector<int> : std::transform(omp_out.begin(), omp_out.end(), omp_in.begin(), omp_out.begin(), std::plus<int>())) initializer(omp_priv = decltype(omp_orig)(omp_orig.size()))

在你使用它的函数之前声明减少，你会很高兴

Answer 2

对于第二个例子，与其存储一个指针然后总是访问同一个元素，不如使用一个局部变量：

    int val = vec[i];
    #pragma omp parallel for reduction (+:val)
    for (int j=0; j<n; ++j)
        val +=j;
    vec[i] = val;

对于第三个循环，我怀疑问题是因为reduction 子句命名了一个变量，但您从未在循环中使用该名称更新该变量，因此编译器不会看到任何减少。使用Eigen 可能会使代码分析起来更加复杂，从而导致循环正常工作。