使用 Eigen::VectorXd 减少 OpenMP

Question

我正在尝试通过减少 OpenMP 来并行化以下循环；

#define EIGEN_DONT_PARALLELIZE
#include <iostream>
#include <cmath>
#include <string>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Eigenvalues>
#include <omp.h>

using namespace Eigen;
using namespace std;

VectorXd integrand(double E)
{
    VectorXd answer(500000);
    double f = 5.*E + 32.*E*E*E*E;
    for (int j = 0; j !=50; j++)
        answer[j] =j*f;
    return answer;
}

int main()
{
    omp_set_num_threads(4);
    double start = 0.;
    double end = 1.;
    int n = 100;
    double h = (end - start)/(2.*n);

    VectorXd result(500000);
    result.fill(0.);
    double E = start;
    result = integrand(E);
    #pragma omp parallel
    { 
    #pragma omp for nowait
    for (int j = 1; j <= n; j++){
        E = start + (2*j - 1.)*h;
        result = result + 4.*integrand(E);
        if (j != n){
            E = start + 2*j*h;
            result = result + 2.*integrand(E);
        }
    }
    }
    for (int i=0; i <50 ; ++i)
        cout<< i+1 << " , "<< result[i] << endl;

    return 0;
}

并行肯定比不并行要快，但是对于所有 4 个线程，结果变化很大。当线程数设置为 1 时，输出正确。 如果有人能帮助我解决这个问题，我将不胜感激......

我正在使用带有编译标志的 clang 编译器；

clang++-3.8 energy_integration.cpp -fopenmp=libiomp5

如果这是失败，那么我将不得不学习实现Boost::thread 或std::thread...

Answer 1

您的代码没有为 OpenMP 定义自定义归约来归约 Eigen 对象。我不确定 clang 是否支持用户定义的缩减（参见OpenMP 4 spec，第 180 页）。如果是这样，您可以声明减少并将reduction(+:result) 添加到#pragma omp for 行。如果没有，你可以自己修改代码，如下：

VectorXd result(500000); // This is the final result, not used by the threads
result.fill(0.);
double E = start;
result = integrand(E);
#pragma omp parallel
{
    // This is a private copy per thread. This resolves race conditions between threads
    VectorXd resultPrivate(500000);
    resultPrivate.fill(0.);
#pragma omp for nowait// reduction(+:result) // Assuming user-defined reductions aren't allowed
    for (int j = 1; j <= n; j++) {
        E = start + (2 * j - 1.)*h;
        resultPrivate = resultPrivate + 4.*integrand(E);
        if (j != n) {
            E = start + 2 * j*h;
            resultPrivate = resultPrivate + 2.*integrand(E);
        }
    }
#pragma omp critical
    {
        // Here we sum the results of each thread one at a time
        result += resultPrivate;
    }
}

您遇到的错误 (in your comment) 似乎是由于大小不匹配造成的。虽然您的代码本身没有一个微不足道的问题，但不要忘记，当 OpenMP 启动每个线程时，它必须为每个线程初始化一个私有 VectorXd。如果没有提供，默认为VectorXd()（大小为零）。使用此对象时，会发生大小不匹配。 omp declare reduction 的“正确”用法将包括初始化程序部分：

#pragma omp declare reduction (+: VectorXd: omp_out=omp_out+omp_in)\
     initializer(omp_priv=VectorXd::Zero(omp_orig.size()))

omp_priv 是私有变量的名称。它由VectorXd::Zero(...) 初始化。使用omp_orig 指定大小。标准 （第 182 页，第 25-27 行）将其定义为：

特殊标识符omp_orig也可以出现在initializer-clause中，它指的是要减少的原始变量的存储。

在我们的例子中（参见下面的完整示例），这是result。所以result.size() 是 500000 并且私有变量被初始化为正确的大小。

#include <iostream>
#include <string>
#include <Eigen/Core>
#include <omp.h>

using namespace Eigen;
using namespace std;

VectorXd integrand(double E)
{
    VectorXd answer(500000);
    double f = 5.*E + 32.*E*E*E*E;
    for (int j = 0; j != 50; j++)   answer[j] = j*f;
    return answer;
}

#pragma omp declare reduction (+: Eigen::VectorXd: omp_out=omp_out+omp_in)\
     initializer(omp_priv=VectorXd::Zero(omp_orig.size()))

int main()
{
    omp_set_num_threads(4);
    double start = 0.;
    double end = 1.;
    int n = 100;
    double h = (end - start) / (2.*n);

    VectorXd result(500000);
    result.fill(0.);
    double E = start;
    result = integrand(E);

#pragma omp parallel for reduction(+:result)
    for (int j = 1; j <= n; j++) {
        E = start + (2 * j - 1.)*h;
        result += (4.*integrand(E)).eval();
        if (j != n) {
            E = start + 2 * j*h;
            result += (2.*integrand(E)).eval();
        }
    }
    for (int i = 0; i < 50; ++i)
        cout << i + 1 << " , " << result[i] << endl;

    return 0;
}