Используйте Eigen Map в сокращении OpenMP

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 не позволяет (или, по крайней мере, не всем реализациям) использовать сокращение для членов класса, но только для переменных.Таким образом, я делаю сокращение для временной матрицы, и у меня есть эта копия в конце 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 Eigen :: CwiseNullaryOp, Eigen :: Matrix>}' в нескалярный тип 'Eigen :: Map, 0, Eigen :: Stride <0, 0>> 'запрошенный инициализатор (omp_priv = Eigen :: MatrixXd :: Zero (omp_orig.rows (), omp_orig.cols ()))

Я получаю, что неявное преобразование из Eigen::MatrixXd в Eigen::Map<Eigen::MatrixXd> не работает в уменьшении omp, но я не знаю, как заставить это работать.

Заранее спасибо

Редактировать 1: Я забыл упомянуть, что яиспользуйте gcc v5.4 на машине с Ubuntu (пробовал и 16.04, и 18.04)

Edit 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, показанный в примере.Использование сокращения таким образом будет довольно дорогим.Если реальный код выполняет что-то вроде фрагмента кода, где значения добавляются на основе вложенных циклов, подобных этому, сокращения можно избежать, разделив работу таким образом, что:

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.