在 OpenMP 缩减中使用 Eigen Map

Question

我想将 Eigen 矩阵与 OpenMP 归约结合使用。

下面是我如何做的一个小例子（它有效）。对象myclass 具有三个属性（一个特征矩阵，两个与其维度对应的整数）和一个成员函数do_something，它使用我定义的总和的 omp 归约，因为特征矩阵不是标准类型。

#include "Eigen/Core"

class myclass {
public:
    Eigen::MatrixXd m_mat;
    int m_n; // number of rows in m_mat
    int m_p; // number of cols in m_mat

    myclass(int n, int p); // constructor

    void do_something(); // omp reduction on `m_mat`
}

myclass::myclass(int n, int p) {
    m_n = n;
    m_p = p;
    m_mat = Eigen::MatrixXd::Zero(m_n,m_p); // init m_mat with null values
}

#pragma omp declare reduction (+: Eigen::MatrixXd: omp_out=omp_out+omp_in)\
    initializer(omp_priv=MatrixXd::Zero(omp_orig.rows(), omp_orig.cols()))

void myclass::do_something() {
    Eigen::MatrixXd tmp = Eigen::MatrixXd::Zero(m_n, m_p); // temporary matrix
#pragma omp parallel for reduction(+:tmp)
    for(int i=0; i<m_n;i++) {
        for(int l=0; l<m_n; l++) {
            for(int j=0; j<m_p; j++) {
                tmp(l,j) += 10;
            }
        }
    }
    m_mat = tmp;
}

问题： OpenMP 不允许（或至少不是所有实现）对类成员使用归约，而只能对变量使用归约。因此，我对一个临时矩阵进行了缩减，并且我在末尾有这个副本m_mat = tmp，我想避免这种情况（因为m_mat 可能是一个大矩阵，我在我的代码中经常使用这种缩减）。

错误修复：我尝试使用 Eigen Map 以使 tmp 对应于存储在 m_mat 中的数据。因此，我将前面代码中的 omp 归约声明和 do_something 成员函数定义替换为：

#pragma omp declare reduction (+: Eigen::Map<Eigen::MatrixXd>: omp_out=omp_out+omp_in)\
    initializer(omp_priv=MatrixXd::Zero(omp_orig.rows(), omp_orig.cols()))

void myclass::do_something() {
    Eigen::Map<Eigen::MatrixXd> tmp = Eigen::Map<Eigen::MatrixXd>(m_mat.data(), m_n, m_p);
#pragma omp parallel for reduction(+:tmp)
    for(int i=0; i<m_n;i++) {
        for(int l=0; l<m_n; l++) {
            for(int j=0; j<m_p; j++) {
                tmp(l,j) += 10;
            }
        }
    }
}

但是，它不再起作用，并且在编译时出现以下错误：

错误：从'const ConstantReturnType {aka const 特征::CwiseNullaryOp, Eigen::Matrix >}' 到非标量类型 ‘Eigen::Map, 0, Eigen::Stride >’ 请求 初始化器(omp_priv=Eigen::MatrixXd::Zero(omp_orig.rows(), omp_orig.cols()))

我知道从 Eigen::MatrixXd 到 Eigen::Map<Eigen::MatrixXd> 的隐式转换在 omp 减少中不起作用，但我不知道如何使它起作用。

提前致谢

编辑 1：我忘了提到我在 Ubuntu 机器上使用 gcc v5.4（尝试了 16.04 和 18.04）

编辑 2： 我修改了我的示例，因为第一个没有减少。这个例子并不完全是我在我的代码中所做的，它只是一个最小的“愚蠢”例子。

Answer 1

正如@ggael 在他们的回答中提到的，Eigen::Map 不能用于此，因为它需要映射到现有存储。如果你确实让它工作，所有线程都将使用相同的底层内存，这会产生竞争条件。

避免您在初始线程中创建的临时变量的最可能解决方案是将成员变量绑定到引用，该引用应该始终有效地用于归约。看起来像这样：

void myclass::do_something() {
    Eigen::MatrixXd &loc_ref = m_mat; // local binding
#pragma omp parallel for reduction(+:loc_ref)
    for(int i=0; i<m_n;i++) {
        for(int l=0; l<m_n; l++) {
            for(int j=0; j<m_p; j++) {
                loc_ref(l,j) += 10;
            }
        }
    }
    // m_mat = tmp; no longer necessary, reducing into the original
}

也就是说，请注意，这仍然会在每个线程中创建零矩阵的本地副本，就像示例中显示的 @ggael 一样。以这种方式使用减少将非常昂贵。如果实际代码正在执行类似于代码 sn-p 的操作，其中值是基于这样的嵌套循环添加的，则可以通过划分工作来避免减少：

1. 每个线程接触矩阵的不同部分
2. 原子用于更新单个值

解决方案 1 示例：

void myclass::do_something() {
    // loop transposed so threads split across l
#pragma omp parallel for
    for(int l=0; l<m_n; l++) {
        for(int i=0; i<m_n;i++) {
            for(int j=0; j<m_p; j++) {
                loc_ref(l,j) += 10;
            }
        }
    }
}

解决方案 2 示例：

void myclass::do_something() {
#pragma omp parallel for
    for(int i=0; i<m_n;i++) {
        for(int l=0; l<m_n; l++) {
            for(int j=0; j<m_p; j++) {
                auto &target = m_mat(l,j);
                // use the ref to get a variable binding through the operator()
                #pragma omp atomic
                target += 10;
            }
        }
    }
}

Answer 2

问题是Eigen::Map 只能在现有内存缓冲区上创建。在您的示例中，底层 OpenMP 实现将尝试执行以下操作：

Eigen::Map<MatrixXd> tmp_0 = MatrixXd::Zero(r,c);
Eigen::Map<MatrixXd> tmp_1 = MatrixXd::Zero(r,c);
...
/* parallel code, thread #i accumulate in tmp_i */
...
tmp = tmp_0 + tmp_1 + ...;

像Map<MatrixXd> tmp_0 = MatrixXd::Zero(r,c) 这样的东西当然是不可能的。 omp_priv 必须是 MatrixXd。我不知道是否可以自定义 OpenMP 创建的私有临时对象的类型。如果不是，您可以通过创建std::vector<MatrixXd> tmps[omp_num_threads()]; 并自己进行最终缩减来手动完成这项工作，或者更好：不要费心制作一个额外的副本，与所有其他工作和副本相比，它在很大程度上可以忽略不计OpenMP 本身。