[Rcpp-devel] Parallel R::bessel_k function in Rcpp

[Hoang Nguyen]

Hi,
I am working with Rcpp and need to call function bessel_k from *R::bessel_k*
in parallel. I got usually crash in Rstudio, and still not find out which
is the mistake

Here is my code in R:

Sys.setenv("PKG_CXXFLAGS"="-fopenmp")
Sys.setenv("PKG_LIBS"="-fopenmp")
library(Rcpp)
sourceCpp("Bessel.cpp")
xseq = rt(10000,df = 8)
param <- c(0.0000000,  2.8284271, -0.4954229,  8.0000000)
BasselCm <- BasselC(xseq,param)

here is my code in Bessel.cpp:
if I set num threads(2) it was crash but run smooth if I set back to 1.

Thank you so much for your consideration.

Best regards.




> sessionInfo()
R version 3.2.4 Revised (2016-03-16 r70336)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
LC_TIME=es_ES.UTF-8
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8
LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=es_ES.UTF-8       LC_NAME=C
LC_ADDRESS=C
[10] LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8
LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] Rcpp_0.12.3

loaded via a namespace (and not attached):
[1] tools_3.2.4               RcppArmadillo_0.6.600.4.0​


--

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

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

using namespace arma ;
using namespace Rcpp ;

// [[Rcpp::export]]
List BasselC(SEXP dataxseq_, SEXP param_) {
  BEGIN_RCPP

  int seed = 126;
  std::srand(seed);

  // Init hyperparams
  Rcpp::NumericVector param(param_);
  double mu = param[0];
  double delta = param[1];
  double beta = param[2];
  double nu = param[3];
  // Rcout <<  "mu " <<  mu <<  " delta " <<  delta <<  " beta " <<  beta <<  " nu " <<  nu  << std::endl;  
  
  
  arma::vec xseq = as<arma::vec>(dataxseq_);
  arma::vec Besselxseq = xseq;
  
  int maxrows = xseq.n_rows;
  
    #ifdef _OPENMP
      omp_set_num_threads(2);
    #endif
      
  #pragma omp parallel for default(none) shared(xseq,mu,delta,beta,nu,maxrows,Besselxseq)    
  for (int i =0; i < maxrows; i++){
    for (int j =0; j < 5000; j++){           //Only for making more computation
      Besselxseq(i) =  log(R::bessel_k(sqrt(pow(beta,2)*(pow(delta,2) + pow(xseq(i) - mu,2))),
                                          (nu + 1)/2, 2)) ;
    }
  }
                      
  
  List z            = List::create(Rcpp::Named("Besselxseq") = Besselxseq);
  return z;
  PutRNGstate();

  END_RCPP
}
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel