Rcpp вызывает segfault RcppArmadillo не

Question

В настоящее время я пытаюсь распараллелить существующую иерархическую схему выборки MCMC. Большая часть моего (теперь последовательного) исходного кода написана в RcppArmadillo, поэтому я хотел бы придерживаться этой структуры для распараллеливания.

Прежде чем начать с распараллеливания моего кода, я прочитал пару сообщений в блоге на Rcpp/Openmp. В большинстве этих сообщений в блоге (например, Drew Schmidt, wrathematics) авторы предупреждают о проблеме безопасности потоков, структур данных R/Rcpp и Openmp. В нижней строке всех сообщений, которые я прочитал до сих пор, R и Rcpp не являются потокобезопасными, не вызывайте их из параллельной прагмы omp.

Из-за этого, в следующем примере Rcpp вызывает Segfault, когда вызывается из R:

#include <Rcpp.h> 
#include <omp.h> 

using namespace Rcpp; 

double rcpp_rootsum_j(Rcpp::NumericVector x) 
{ 
    Rcpp::NumericVector ret = sqrt(x); 
    return sum(ret); 
} 

// [[Rcpp::export]] 
Rcpp::NumericVector rcpp_rootsum(Rcpp::NumericMatrix x, int cores = 2) 
{ 
    omp_set_num_threads(cores); 
    const int nr = x.nrow(); 
    const int nc = x.ncol(); 
    Rcpp::NumericVector ret(nc); 

    #pragma omp parallel for shared(x, ret) 
    for (int j=0; j<nc; j++) 
    ret[j] = rcpp_rootsum_j(x.column(j)); 

    return ret; 
} 

Как Дрю объясняет в своем блоге, то выдаёт ошибку сегментации происходит из-за «скрытой» копии, что делает Rcpp в обращении к ret[j] = rcpp_rootsum_j(x.column(j));.

Поскольку я заинтересован в поведении RcppArmadillo в случае распараллеливания, у меня есть пример преобразованного Дрю:

//[[Rcpp::depends(RcppArmadillo)]] 
#include <RcppArmadillo.h> 
#include <omp.h> 

double rcpp_rootsum_j_arma(arma::vec x) 
{ 
    arma::vec ret = arma::sqrt(x); 
    return arma::accu(ret); 
} 

// [[Rcpp::export]] 
arma::vec rcpp_rootsum_arma(arma::mat x, int cores = 2) 
{ 
    omp_set_num_threads(cores); 
    const int nr = x.n_rows; 
    const int nc = x.n_cols; 
    arma::vec ret(nc); 

    #pragma omp parallel for shared(x, ret) 
    for (int j=0; j<nc; j++) 
    ret(j) = rcpp_rootsum_j_arma(x.col(j)); 

    return ret; 
} 

Интересно, что семантически эквивалентный код не вызывает Segfault.

Вторая вещь, которую я заметил во время моего исследования в том, что вышеупомянутое заявление (R и Rcpp не поточно, не называйте их изнутри OMP параллельно прагме), кажется, не всегда имеет место, чтобы быть правдой , Например, вызов в следующем примере не вызывает segfault, хотя мы читаем и записываем структуры данных Rcpp.

#include <Rcpp.h> 
#include <omp.h> 

// [[Rcpp::export]] 
Rcpp::NumericMatrix rcpp_sweep_(Rcpp::NumericMatrix x, Rcpp::NumericVector vec) 
{ 
    Rcpp::NumericMatrix ret(x.nrow(), x.ncol()); 

    #pragma omp parallel for default(shared) 
    for (int j=0; j<x.ncol(); j++) 
    { 
    #pragma omp simd 
    for (int i=0; i<x.nrow(); i++) 
     ret(i, j) = x(i, j) - vec(i); 
    } 

    return ret; 
} 

Мои вопросы

1. Почему код из первого примера вызывает Segfault, но код из примера два и три не?
2. Как узнать, может ли безопасный вызов метода (arma::mat.col(i)) или если небезопасно вызывать метод (Rcpp::NumericMatrix.column(i))? Должен ли я каждый раз читать исходный код среды?
3. Любые предложения о том, как избежать таких «непрозрачных» ситуаций, как в примере один?

Это может быть чистое совпадение, что мой пример RcppArmadillo не подведет. См. Комментарии Dirks ниже.

EDIT 1

В своем ответе и в обоих своих комментариях Dirk настоятельно рекомендует более внимательно изучить примеры в Rcpp галерее.

Вот мои первоначальные предположения:

1. Извлечение строк, столбцов и т.д. внутри прагма OpenMP, как правило, не поточно, так как он может перезвонить в R, чтобы выделить новое пространство в памяти для скрытой копии.
2. Поскольку RcppArmadillo полагается на ту же модель облегченного/прокси для структур данных, что и Rcpp, мое первое предположение также относится к RcppArmadillo.
3. Структуры данных из пространства имен std должны быть как можно более безопасными, поскольку они не используют одну и ту же облегченную/прокси-схему.
4. Примитивные типы данных также не должны вызывать проблем, поскольку они живут в стеке и не требуют, чтобы R выделяла и управляла памятью.

Optimizing Code vs...

arma::mat temp_row_sub = temp_mat.rows(x-2, x+2); 

Hierarchical Risk Parity...

interMatrix(_, i) = MAT_COV(_, index_asset); // 3rd code example 3rd method 

Using RcppProgress...

thread_sum += R::dlnorm(i+j, 0.0, 1.0, 0); // subsection OpenMP support 

На мой взгляд, первый и второй пример явно мешают мои предположения, сделанные в точке о ne и два. Пример три также дает мне головную боль, так как для меня это выглядит как вызов к R ...

Мои обновленные вопросы

1. Где разница между, например один/два, и мой первый пример кода?
2. Где я заблудился в своих предположениях?

Любые рекомендации, помимо RcppGallery и GitHub, о том, как лучше понять взаимодействие Rcpp и OpenMP?

Answer 1

Прежде чем начать с распараллеливания кода, я прочитал пару сообщений в блоге на Rcpp/Openmp. В большинстве этих сообщений в блоге (например, Дрю Шмидт, wrathematics) авторы предупреждают о проблеме безопасности потоков, структур данных R/Rcpp и Openmp. В нижней строке всех сообщений, которые я прочитал до сих пор, R и Rcpp не являются потокобезопасными, не вызывайте их из параллельной прагмы omp.

Это известное ограничение самого R, не являющегося потокобезопасным. Это означает, что вы не можете перезвонить или запустить R-события, что может произойти с Rcpp, если вы не будете осторожны. Чтобы быть более простым: ограничение не имеет ничего общего с Rcpp, это просто означает, что вы не можете слепо заглянуть в OpenMP через Rcpp. Но вы можете, если будете осторожны.

Мы бесчисленных примеров успеха с помощью OpenMP и связанных с ними инструментами, как в многочисленных пакетах на CRAN, на Rcpp галерее и через расширение пакеты, такие как RcppParallel.

Вы, кажется, были очень избирательно в том, что вы выбрали для чтения на эту тему, и в итоге у вас что-то произошло между неправильным и вводящим в заблуждение. Я предлагаю вам обратиться к нескольким примерам на Rcpp Gallery, которые касаются OpenMP/RcppParallel, поскольку они имеют дело с самой проблемой. Или если вы спешите: посмотрите RVector и RMatrix в документации RcppParallel.

Ресурсы:

• Six posts on OpenMP at RcppGallery
• RcppParallel site
• RcppParallel CRAN package

и ваш самый большой ресурс может быть какой-то целевой поиск в GitHub для кода с участием R, C++ и OpenMP. Это приведет вас к многочисленным рабочим примерам.