使用 Rcpp/Armadillo 和具有用户定义函数的 openMP praralel 时出现段错误

Question

我试图将 rcpp/armadillo 与 openmp 一起使用来加速 R 中的循环。循环采用一个矩阵，其中每一行包含位置向量的索引（如果是二维位置，则为矩阵）作为输入（和其他 matrix/vec 被使用）。在循环中，我提取了输入索引矩阵的每一行并找到相应的位置，计算距离矩阵和协方差矩阵，进行 cholesky 和 backsolve，将 backsolve 结果保存到一个新矩阵。这是 rcpp 代码：

`#include <iostream>
#include <RcppArmadillo.h>
#include <omp.h>
#include <Rcpp.h>

// [[Rcpp::plugins(openmp)]]
using namespace Rcpp;
using namespace arma;
using namespace std;

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
mat NZentries_new2 (int m, int nnp, const mat& locs, const umat& revNNarray, const mat& revCondOnLatent, const vec& nuggets, const vec covparms){
// initialized the output matrix
mat Lentries=zeros(nnp,m+1);
// initialized objects in parallel part
int n0; //number of !is_na elements
uvec inds;//
vec revCon_row;//
uvec inds00;//
vec nug;//
mat covmat;//
vec onevec;//
vec M;//
mat dist;//
int k;//
omp_set_num_threads(2);// selects the number of cores to use.
#pragma omp parallel for shared(locs,revNNarray,revCondOnLatent,nuggets,nnp,m,Lentries) private(k,M,dist,onevec,covmat,nug,n0,inds,revCon_row,inds00) default(none) schedule(static)
for (k = 0; k < nnp; k++) {
// extract a row to work with
inds=revNNarray.row(k).t();
revCon_row=revCondOnLatent.row(k).t();

if (k < m){
  n0=k+1;
} else {
  n0=m+1;
}
// extract locations
inds00=inds(span(m+1-n0,m))-ones<uvec>(n0);

nug=nuggets.elem(inds00) % (ones(n0)-revCon_row(span(m+1-n0,m))); // vec is vec, cannot convert to mat
dist=calcPWD2(locs.rows(inds00));

#pragma omp critical
{
  //calculate covariance matrix
  covmat= MaternFun(dist,covparms) + diagmat(nug) ; // summation from arma
}

// get last row of inverse Cholesky
onevec = zeros(n0);
onevec[n0-1] = 1;
M=solve(chol(covmat,"upper"),onevec);
// save the entries to matrix
Lentries(k,span(0,n0-1)) = M.t();
}
return Lentries;
}`

当前版本可以正常运行，但速度较慢（几乎与没有并行版本相同），如果我将omp critical bracket 中的线去掉，会导致段错误，R 会崩溃。这个 MaterFun 是我定义的一个函数，和其他几个小函数一起定义如下。所以我的问题是为什么 MaternFun 必须留在关键部分。

// [[Rcpp::export]]
    mat MaternFun( mat distmat, vec covparms ){

      int d1 = distmat.n_rows;
      int d2 = distmat.n_cols;
      int j1;
      int j2;
      mat covmat(d1,d2);
      double scaledist;

      double normcon = covparms(0)/(pow(2.0,covparms(2)-1)*Rf_gammafn(covparms(2)));

      for (j1 = 0; j1 < d1; j1++){
        for (j2 = 0; j2 < d2; j2++){
          if ( distmat(j1,j2) == 0 ){
            covmat(j1,j2) = covparms(0);
          } else {
            scaledist = distmat(j1,j2)/covparms(1);
            covmat(j1,j2) = normcon*pow( scaledist, covparms(2) )*
              Rf_bessel_k(scaledist,covparms(2),1.0);
          }
        }
      }
      return covmat;
    }

    // [[Rcpp::export]]
    double dist2(double lat1,double long1,double lat2,double long2) {
      double dist = sqrt(pow(lat1 - lat2, 2) + pow(long1 - long2, 2)) ;
      return (dist) ;
    }
   // [[Rcpp::export]]
    mat calcPWD2( mat x) {//Rcpp::NumericMatrix
      int outrows = x.n_rows ;

      int outcols = x.n_rows ;
      mat out(outrows, outcols) ;
      for (int arow = 0 ; arow < outrows ; arow++) {
        for (int acol = 0 ; acol < outcols ; acol++) {
          out(arow, acol) = dist2(x(arow, 0),x(arow, 1),
                                  x(acol, 0),x(acol, 1)) ; //extract element from mat 
        }
      }
      return (out) ;
    }

这是一些用于测试 R 中的 MaterFun 的示例输入：library(fields) distmat=rdist(1:5) # distance matrix covparms=c(1,0.2,1.5)

Answer 1

问题是有两次调用 R 数学函数（Rf_bessel_k 和 Rf_gammafn）需要单线程而不是并行访问.

为了解决这个问题，让我们在 C++ 中通过 BH to obtain the cyl_bessel_k and tgamma functions. Alternatively, there is always the option of reimplementing R's besselK and gamma 添加对 boost 的依赖，这样它就不会使用单线程 R变体。

这给出：

#include <Rcpp.h>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double besselK_boost(double x, double v) {
  return boost::math::cyl_bessel_k(v, x);
}

// [[Rcpp::export]]
double gamma_fn_boost(double x) {
  return boost::math::tgamma(x);
}

测试代码

x0 = 9.536743e-07
nu = -10
all.equal(besselK(x0, nu), besselK_boost(x0, nu))
# [1] TRUE

x = 2
all.equal(gamma(x), gamma_fn_boost(x))
# [1] TRUE

注意：boost 变体的参数顺序与 R 的不同：

cyl_bessel_k(v, x)
Rf_bessel_k(x, v, expon.scaled = FALSE)

从这里，我们可以修改MaternFun。不幸的是，因为缺少 calcPWD2，我们最多只能切换到使用 boost 并合并到 OpenMP protections.

#include <RcppArmadillo.h>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

#ifdef _OPENMP
#include <omp.h>
#endif

// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::depends(BH)]]
// [[Rcpp::plugins(openmp)]]

// [[Rcpp::export]]
arma::mat MaternFun(arma::mat distmat, arma::vec covparms) {

  int d1 = distmat.n_rows;
  int d2 = distmat.n_cols;
  int j1;
  int j2;
  arma::mat covmat(d1,d2);
  double scaledist;

  double normcon = covparms(0) /
    (pow(2.0, covparms(2) - 1) * boost::math::tgamma(covparms(2)));

  for (j1 = 0; j1 < d1; ++j1){
    for (j2 = 0; j2 < d2; ++j2){
      if ( distmat(j1, j2) == 0 ){
        covmat(j1, j2) = covparms(0);
      } else {
        scaledist = distmat(j1, j2)/covparms(1);
        covmat(j1, j2) = normcon * pow( scaledist, covparms(2) ) *
          boost::math::cyl_bessel_k(covparms(2), scaledist);
      }
    }
  }
  return covmat;
}

使用 Rcpp/Armadillo 和具有用户定义函数的 openMP praralel 时出现段错误

Segment fault when using Rcpp/Armadillo and openMP prarallel with user-defined function

r

openmp

rcpp