Re: [eigen] Use Eigen Map in OpenMP reduction

[Gael Guennebaud <gael.guennebaud@xxxxxxxxx>]

Hi,

I answered directly on SO as this will make my explanation more visible to future readers.

cheers,

Gael.

On Thu, Jan 10, 2019 at 7:30 PM 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