Rcpp 导致段错误 RcppArmadillo 不会

Question

我目前正在尝试并行化现有的分层 MCMC 采样方案。我的大部分（现在是连续的）源代码都是用 RcppArmadillo 编写的，所以我也想坚持使用这个并行化框架。

在开始并行化我的代码之前，我阅读了 Rcpp/Openmp 上的一些博客 post。在这些博客的大多数 post 中（例如 Drew Schmidt, wrathematics），作者警告线程安全、R/Rcpp 数据结构和 Openmp 问题。到目前为止我读过的所有 posts 的底线是，R 和 Rcpp 不是线程安全的，不要从 omp parallel pragma 中调用它们。

因此，当从 R 调用时，以下 Rcpp 示例会导致段错误：

#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 在他的博客 post 中解释的那样，段错误的发生是由于 "hidden" 副本，Rcpp 在调用 ret[j] = rcpp_rootsum_j(x.column(j));.

时产生的

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

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

有趣的是，语义上等效的代码不会导致段错误。

我在研究过程中注意到的第二件事是，上述声明（R 和 Rcpp 不是线程安全的，不要从 omp parallel pragma 中调用它们) 似乎并不总是成立。例如，下一个示例中的调用不会导致段错误，尽管我们正在读取和写入 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;
}

我的问题

~~为什么第一个示例中的代码会导致段错误，而示例二和三中的代码却不会？~~

~~我怎么知道调用方法是安全的 (arma::mat.col(i)) 还是调用方法不安全 (Rcpp::NumericMatrix.column(i))？每次都要看框架的源码吗？~~

~~关于如何避免这些 "opaque" 情况的任何建议，如示例 1？~~

我的 RcppArmadillo 示例没有失败可能纯属巧合。请参阅下面的 Dirks 评论。

编辑 1

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

以下是我的初步假设：

在 OpenMp pragma 中提取行、列等通常不是线程安全的，因为它可能会回调 R 以在内存中为隐藏副本分配新的 space。
因为 RcppArmadillo 依赖与 Rcpp 相同的 lightweight/proxy 数据结构模型，所以我的第一个假设也适用于 RcppArmadillo。
std namespace 中的数据结构在某种程度上应该更安全，因为它们不使用相同的 lightweight/proxy 方案。
原始数据类型也不应该引起问题，因为它们存在于堆栈中并且不需要 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

在我看来，第一个和第二个示例明显干扰了我在第一点和第二点中所做的假设。示例三也让我头疼，因为对我来说它看起来像是对 R 的调用...

我更新的问题

例子 one/two 和我的第一个代码例子有什么区别？
我在哪里迷失了自己的假设？

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

Answer 1

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。但如果你小心的话，你可以。

我们有无数个成功 OpenMP 和相关工具的例子，在 CRAN、Rcpp Gallery 和扩展中的许多包中像 RcppParallel 这样的软件包。

您似乎非常有选择性地选择阅读有关该主题的内容，并且最终得出的结果介于错误和误导之间。我建议您转向 Rcpp Gallery 上处理 OpenMP / RcppParallel 的几个示例，因为它们处理的正是这个问题。或者，如果您赶时间：在 RcppParallel 文档中查找 RVector 和 RMatrix。

资源：

您最大的资源可能是在 GitHub 上进行一些有针对性的搜索，以查找涉及 R、C++ 和 OpenMP 的代码。它将引导您找到许多工作示例。

Rcpp 导致段错误 RcppArmadillo 不会

Rcpp causes segfault RcppArmadillo does not

parallel-processing

openmp

armadillo

rcpp