parallel-processing - Rcpp 导致段错误 RcppArmadillo 不会

Question

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

在开始并行化我的代码之前，我阅读了几篇关于 Rcpp/Openmp 的博客文章。在这些博客文章的大部分(例如 Drew Schmidt, wrathematics)中，作者警告线程安全、R/Rcpp 数据结构和 Openmp 问题。到目前为止我读过的所有帖子的底线是，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 在他的博客文章中所解释的那样，段错误的发生是由于“隐藏”副本，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;
}

我的问题

1. 为什么第一个示例中的代码会导致段错误，而示例二和三中的代码却不会？
<罢工>
2. 我怎么知道调用方法是安全的 (arma::mat.col(i)) 还是调用方法不安全 (Rcpp::NumericMatrix.column(i))？我是否必须每次都阅读框架的源代码？
3. 关于如何避免示例一中的这些“不透明”情况的任何建议？

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

编辑 1

在他的回答和他的两条评论中，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

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

我们在 CRAN、Rcpp Gallery 和扩展包(如 RcppParallel)上的许多包中都有 无数成功的 OpenMP 和相关工具示例。

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

资源:

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

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