Rcpp 会导致段错误 RcppArmadillo 不会

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

在开始并行化我的代码之前，我阅读了几篇关于 Rcpp/Openmp 的博客文章。在大多数这些博客文章中（例如德鲁·施密特，wrathematics）作者对线程安全、R/Rcpp 数据结构和 Openmp 问题发出警告。到目前为止我读过的所有帖子的底线是，R 和 Rcpp 不是线程安全的，不要从 omp 并行编译指示中调用它们。

因此，当从 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 在他的博客文章中解释的那样，段错误是由于 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 并行编译指示中调用它们）似乎并不总是成立。例如，下一个示例中的调用不会导致段错误，尽管我们正在读取和写入 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. 为什么第一个示例中的代码会导致段错误，而示例二和示例三中的代码却不会？ 罢工>
2. 我怎么知道调用方法是否安全（arma::mat.col(i)）或者如果调用方法不安全（Rcpp::NumericMatrix.column(i)）？我每次都必须阅读框架的源代码吗？
3. 关于如何避免这些“不透明”情况（如示例一）的任何建议？

在开始并行化我的代码之前，我阅读了几篇关于 Rcpp/Openmp 的博客文章。在大多数这些博客文章（例如 Drew Schmidt、wrathematics）中，作者对线程安全、R/Rcpp 数据结构和 Openmp 问题提出警告。到目前为止我读过的所有帖子的底线是，R 和 Rcpp 不是线程安全的，不要从 omp 并行编译指示中调用它们。

这是 R 本身不是线程安全的一个众所周知的限制。这意味着您无法回调或触发 R 事件——除非您小心，否则 Rcpp 可能会发生这种情况。更简单地说：该约束与 Rcpp 无关，它只是意味着您不能通过 Rcpp 盲目地进入 OpenMP。但如果你小心的话，你可以做到。

We have 无数示例successOpenMP 和相关工具均位于 CRAN、Rcpp Gallery 上的众多软件包中以及通过 RcppParallel 等扩展包。

你似乎已经非常有选择性在你选择阅读的关于这个主题的内容中，你最终得到了一些介于错误和误导之间的东西。我建议你看看上面的几个例子RCPP画廊它处理 OpenMP / RcppParallel 因为它们处理这个问题。或者，如果您赶时间：请向上查找RVector and RMatrix在 RcppParallel 文档中。

资源：

• RcppGallery 上有 6 个关于 OpenMP 的帖子
• Rcpp并行站点
• Rcpp并行CRAN包

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