Re: [eigen] Use Eigen Map in OpenMP reduction

[Francois Fayard <fayard@xxxxxxxxxxxxx>]

It needs to be confirmed by a close look at the OpenMP standard, but I don’t think it allows to reduce into anything else than a plain old data. Or at least, I don’t think that any implementation of OpenMP allows that.

That’s not answering the question, but TBB allows to define your own reduction and what you want to do can be achieved with it.

François

On 10 Jan 2019, at 17:12, Ghislain Durif <gd.dev@xxxxxxxxxxxxxxx> wrote:

Hi,

I have a question regarding the use of Eigen matrices combined with OpenMP reduction. I posted the same question on Stackoverflow (https://stackoverflow.com/questions/53957734/use-eigen-map-in-openmp-reduction) with no clear answer so far. I copy here my entire question with a little "dumb" example. I hope someone will be able to help me.

How to use Eigen matrices in combination with OpenMP reduction ? (it is complicated to do something efficient otherwise in the algorithm that I implement)

In the following is a small example of how I do it (and it works). The object myclass has three attributes (an Eigen matrix, two integers corresponding to its dimension) and a member function do_something that uses an omp reduction on a sum which I define because Eigen matrices are not standard types.

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

Problem: OpenMP does not allow (or at least not all implementations) to use reduction on class members but only on variables. Thus, I do the reduction on a temporary matrix and I have this copy at the end m_mat = tmp which I would like to avoid (because m_mat can be a big matrix and I use this reduction a lot in my code).

Wrong fix: I tried to use Eigen Map so that tmp corresponds to data stored in m_mat. Thus, I replaced the omp reduction declaration and the do_something member function definition in the previous code with:

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

However, it does not work anymore and I get the following error at compilation:

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()))

I get that the implicit conversion from Eigen::MatrixXd to Eigen::Map<Eigen::MatrixXd> does not work in the omp reduction but I don't know how to make it work.

I use gcc v5.4 on a Ubuntu machine (tried both 16.04 and 18.04)


Thanks in advance,
Best regards and happy new year

-- 
Ghislain Durif