使用 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 定义自定义缩减来缩减特征对象。我不确定 clang 是否支持用户定义的缩减（请参阅OpenMP 4 规范，第 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;
    }
}

您收到的错误（在您的评论中）似乎是由于尺寸不匹配造成的。虽然您的代码本身并不重要，但不要忘记，当 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也可以出现在初始化子句中，它将引用要减少的原始变量的存储。

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

#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;
}