在 OpenMP 缩减中使用 Eigen Map

Question

我想结合使用特征矩阵和 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 不允许（或至少不是所有实现）对 class 成员使用归约，但只能对变量使用归约。因此，我对一个临时矩阵进行了缩减，我在末尾有这个副本 m_mat = tmp 我想避免（因为 m_mat 可能是一个大矩阵，我在我的代码）。

错误修复： 我尝试使用 Eigen Map 以便 tmp 对应于存储在 m_mat 中的数据。因此，我将之前代码中的 omp reduction 声明和 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;
            }
        }
    }
}

但是，它不再起作用，我在编译时收到以下错误：

error: conversion from ‘const ConstantReturnType {aka const Eigen::CwiseNullaryOp, Eigen::Matrix >}’ to non-scalar type ‘Eigen::Map, 0, Eigen::Stride<0, 0> >’ requested initializer(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： 我修改了示例，因为第一个示例没有减少。这个例子并不完全是我在我的代码中所做的，它只是一个最小的 "dumb" 例子。

Answer 1

问题是 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 本身完成。

Answer 2

正如@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 一样。以这种方式使用减少将非常昂贵。如果实际代码正在执行类似代码片段的操作，其中基于嵌套循环添加值，则可以通过划分工作来避免减少：

每个线程接触矩阵的不同部分
一个原子用于更新个体值

解决方案 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;
            }
        }
    }
}

在 OpenMP 缩减中使用 Eigen Map

Use Eigen Map in OpenMP reduction

c++

openmp

eigen