Rcpp 导致段错误 RcppArmadillo 没有

Question

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

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

因此，以下 Rcpp 示例在从 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)); 时生成的“隐藏”副本而发生的。

由于我对 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. 关于如何避免这些“不透明”情况（如示例一）有什么建议吗？

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

编辑 1

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

这是我最初的假设：

1. 在 OpenMp pragma 中提取行、列等通常不是线程安全的，因为它可能会回调到 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

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

我更新的问题

1. 示例一/二与我的第一个代码示例的区别在哪里？
2. 我在哪里迷失了自己的假设？

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

Answer 1

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

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

我们有无数个使用 OpenMP 和相关工具的成功示例，包括在 CRAN、Rcpp Gallery 上的众多包中以及通过 RcppParallel 等扩展包。

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

资源：

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

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