Rcpp导致段错RcppArmadillo没有
我目前正在尝试并行化现有的分层MCMC采样方案。我的大部分（现在顺序）源代码都是用RcppArmadillo编写的，所以我也想坚持使用这个框架进行并行化。
在开始并行化我的代码之前，我已经阅读了几篇关于Rcpp / Openmp的博客文章。在大多数博客文章中（例如Drew Schmidt，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;
}
我的问题
为什么第一个示例中的代码会导致段错误，但示例二和三的代码不会？如何调用方法（arma::mat.col(i)）或调用方法（）是否安全？我怎么知道Rcpp::NumericMatrix.column(i)？我每次都必须阅读框架的源代码吗？有关如何避免像示例一样的“不透明”情况的任何建议？
我的RcppArmadillo示例并没有失败，这可能是纯粹的巧合。请参阅下面的Dirks评论。
编辑1
在他的回答和他的两篇评论中，德克强烈建议更仔细地研究Rcpp画廊中的例子。
以下是我最初的假设：
- 在OpenMp pragma中提取行，列等通常不是线程安全的，因为它可能会回调到R以在内存中为隐藏副本分配新空间。
- 因为RcppArmadillo依赖于与Rcpp相同的轻量/代理模型用于数据结构，所以我的第一个假设也适用于RcppArmadillo。
- 来自std命名空间的数据结构应该更安全，因为它们不使用相同的轻量级/代理方案。
- 原始数据类型也不应该引起问题，因为它们存在于堆栈中并且不需要R来分配和管理内存。
arma::mat temp_row_sub = temp_mat.rows(x-2, x+2);
interMatrix(_, i) = MAT_COV(_, index_asset); // 3rd code example 3rd method
thread_sum += R::dlnorm(i+j, 0.0, 1.0, 0); // subsection OpenMP support
在我看来，第一个和第二个例子明显干扰了我在第一点和第二点所作的假设。示例三也让我头疼，因为对我而言，它看起来像是对R的调用......
我的最新问题
- 示例一/二和我的第一个代码示例之间的区别在哪里？
- 我在假设中迷失了哪里？
除了RcppGallery和GitHub之外，还有什么建议可以更好地了解Rcpp和OpenMP的交互？
Stack Overflow用户
发布于 2019-06-03 13:43:53
在开始并行化我的代码之前，我已经阅读了几篇关于Rcpp / Openmp的博客文章。在大多数博客文章中（例如Drew Schmidt，wrathematics），作者警告线程安全，R / Rcpp数据结构和Openmp问题。我到目前为止读过的所有帖子的底线是，R和Rcpp不是线程安全的，不要在omp并行编译指示中调用它们。
这是R本身不为线程安全的众所周知的限制。这意味着你不能回电，或触发R事件 - 除非你小心，否则可能会发生Rcpp事件。更简单：约束与Rcpp无关，它只是意味着你不能盲目地通过Rcpp进入OpenMP。但是如果你小心的话，你可以。
我们在CRAN的众多软件包，Rcpp Gallery和RcppParallel等扩展软件包中都有无数的OpenMP和相关工具成功案例。
你似乎对你在这个主题上所选择的内容非常挑剔，而你最终得到的东西介于错误和误导之间。我建议你转向Rcpp Gallery上的几个例子来处理OpenMP / RcppParallel，因为它们处理了这个问题。或者如果你赶时间：查看RVector并RMatrix在RcppParallel文档中。
资源：
并且您最大的资源可能是在GitHub上针对涉及R，C ++和OpenMP的代码进行一些有针对性的搜索。它将引导您进行大量工作示例。
https://stackoverflow.com/questions/-100006893
复制相似问题