Rcpp导致段错误RcppArmadillo不会 [英] Rcpp causes segfault RcppArmadillo does not

问题描述

我目前正在尝试并行化现有的分层MCMC采样方案.我的大部分(现在是顺序的)源代码都是用RcppArmadillo编写的，因此，我也想坚持使用此框架进行并行化.

I'm currently trying to parallelize an existing hierarchical MCMC sampling scheme. The majority of my (by now sequential) source code is written in RcppArmadillo, so I'd like to stick with this framework for parallelization, too.

在开始并行化代码之前，我已经阅读了Rcpp/Openmp上的几篇博客文章.在大多数这些博客文章中(例如， wrewthematics的德鲁·施密特(arew Schmidt))，作者都对以下问题提出了警告.线程安全，R/Rcpp数据结构和Openmp.到目前为止，我读过的所有帖子的底线是R和Rcpp都不是线程安全的，请不要从omp并行编译指示中调用它们.

Before starting with parallelizing my code I have read a couple of blog posts on Rcpp/Openmp. In the majority of these blog posts (e.g. Drew Schmidt, wrathematics) the authors warn about the issue of thread safety, R/Rcpp data structures and Openmp. The bottom line of all posts I have read so far is, R and Rcpp are not thread safe, don't call them from within an omp parallel pragma.

因此，下面的Rcpp示例从R调用时会导致段错误:

Because of this, the following Rcpp example causes a segfault, when called from 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;
}

如Drew在他的博客文章中所述，该段错误是由于Rcpp在调用ret[j] = rcpp_rootsum_j(x.column(j));时创建的隐藏"副本而发生的.

As Drew explains in his blog post, the segfault happens due to a "hidden" copy, which Rcpp makes in the call to ret[j] = rcpp_rootsum_j(x.column(j));.

由于我对并行化情况下的RcppArmadillo的行为感兴趣，因此我转换了Drew的示例:

Since I'm interested in the behavior of RcppArmadillo in case of parallelization, I have converted Drew's example:

//[[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;
}

有趣的是，语义上等效的代码不会引起段错误.

Interestingly the semantically equivalent code does not causes a segfault.

我在研究过程中注意到的第二件事是，上述陈述( R和Rcpp都不是线程安全的，不要从omp并行编译指示内调用它们)似乎并不总是坚持是真实的.例如，尽管我们正在读取和写入Rcpp数据结构，但下一个示例中的调用不会导致段错误.

The second thing I have noticed during my research is, that the aforementioned statement (R and Rcpp are not thread safe, don't call them from within an omp parallel pragma) seems to not always hold to be true. For example the call in the next example does not cause a segfault, although we're reading and writing to Rcpp data structures.

#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. 为什么第一个示例中的代码会导致段错误，而第二个示例和三个示例中的代码却不会?
2. 如何知道调用方法(arma::mat.col(i))的安全性还是调用方法(Rcpp::NumericMatrix.column(i))的安全性?我是否每次都必须阅读框架的源代码?
3. 关于如何避免出现示例一中的不透明"情况的任何建议?

1. Why does the code from the first example causes a segfault but the code from example two and three not?
2. How do I know, if it is safe to call a method (arma::mat.col(i)) or if it is unsafe to call a method (Rcpp::NumericMatrix.column(i))? Do I have to read the framework's source code every time?
3. Any suggestions on how to avoid these "opaque" situations like in example one?

我的RcppArmadillo示例不会失败可能纯粹是巧合.请参阅下面的Dirks评论.

编辑1

在他的回答和他的两条评论中，Dirk强烈建议更仔细地研究Rcpp Gallery中的示例.

In his answer and in both of his comments Dirk strongly recommends to study more closely the examples in the Rcpp Gallery.

这是我最初的假设:

1. 在OpenMp编译指示中提取行，列等通常不是线程安全的，因为它可能会回调R来为隐藏副本分配内存中的新空间.
2. 由于RcppArmadillo依赖与Rcpp相同的数据结构轻量级/代理模型，因此我的第一个假设也适用于RcppArmadillo.
3. std名称空间中的数据结构应该以某种方式更安全，因为它们不使用相同的轻量级/代理方案.
4. 原始数据类型也不会引起问题，因为它们存在于堆栈中，并且不需要R来分配和管理内存.

优化代码与...

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

分级风险均等...

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

使用RcppProgress ...

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

在我看来，第一个和第二个例子显然干扰了我在第一点和第二点所作的假设.示例三也让我头疼，因为对我而言，这似乎是对R ...的调用.

In my opinion the first and second example clearly interfere with my assumptions made in point one and two. Example three also gives me headaches, since for me it looks like a call to R...

我更新的问题

1. 示例一/二和我的第一个代码示例有何区别?
2. 我的假设在哪里迷路了?

除RcppGallery和GitHub以外，还有关于如何更好地了解Rcpp和OpenMP交互的任何建议吗?

Any recommendations, besides RcppGallery and GitHub, on how to get a better idea of the interaction of Rcpp and OpenMP?

推荐答案

在开始并行化代码之前，我已经阅读了Rcpp/Openmp上的几篇博客文章.在大多数这些博客帖子中(例如，wrewthematics的Drew Schmidt)，作者警告线程安全，R/Rcpp数据结构和Openmp的问题.到目前为止，我读过的所有帖子的底线是R和Rcpp都不是线程安全的，请不要从omp并行编译指示中调用它们.

Before starting with parallelizing my code I have read a couple of blog posts on Rcpp/Openmp. In the majority of these blog posts (e.g. Drew Schmidt, wrathematics) the authors warn about the issue of thread safety, R/Rcpp data structures and Openmp. The bottom line of all posts I have read so far is, R and Rcpp are not thread safe, don't call them from within an omp parallel pragma.

这是R本身不是线程安全的众所周知的限制.这意味着您不能回叫或触发R事件-除非您小心，否则Rcpp可能会发生.更简单地说:约束与Rcpp无关，它只是意味着您不能通过Rcpp盲目地进入OpenMP.但是如果您小心的话就可以.

That is a well-known limitation of R itself not being thread-safe. That means you cannot call back, or trigger R events -- which may happen with Rcpp unless you are careful. To be more plain: The constraint has nothing to do with Rcpp, it simply means you cannot blindly drop into OpenMP via Rcpp. But you can if you're careful.

在CRAN上的许多程序包，Rcpp Gallery上以及通过扩展程序包(例如RcppParallel)，我们都有带有OpenMP和相关工具的成功的无数示例.

We have countless examples of success with OpenMP and related tools both in numerous packages on CRAN, on the Rcpp Gallery and via extension packages like RcppParallel.

在选择阅读该主题时，您似乎非常有选择性，结果导致错误和误导之间存在某种差异.我建议您转到 Rcpp Gallery 上的几个示例，它们处理OpenMP/RcppParallel，因为它们处理了很有问题.或者，如果您需要赶时间:在RcppParallel文档中查找RVector和RMatrix.

You appear to have been very selective in what you chose to read on this topic, and you ended up with something somewhere between wrong and misleading. I suggest you turn to the several examples on the Rcpp Gallery which deal with OpenMP / RcppParallel as they deal with the very problem. Or if you're in a hurry: look up RVector and RMatrix in the RcppParallel documentation.

资源:

• 在RcppGallery上OpenMP上的六篇文章
• RcppParallel网站
• RcppParallel CRAN软件包

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

和您最大的资源可能是在GitHub上针对性地搜索涉及R，C ++和OpenMP的代码.它将带您找到许多可行的示例.

and your largest resource may be some targeted search at GitHub for code involving R, C++ and OpenMP. It will lead you to numerous working examples.

这篇关于Rcpp导致段错误RcppArmadillo不会的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！