我目前正在尝试并行化现有的分层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画廊中的例子。
以下是我最初的假设:
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的交互?
发布于 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
复制相似问题