CODEWITHSUNDEEP
rcpp
Ignacio
Ignacio

Reputation: 7928

Rstudio crashes with Rcpp and OpenMP function

This is a follow up question to dqrng with Rcpp for drawing from a normal and a binomial distribution. I tried to implement the answer but instead of drawing from a single distribution I'm drawing from 3. This is the code that I wrote:

// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
  dqrng::xoshiro256plus rng(42);
  arma::mat out(n*m,3);
  // ok to use rng here
#pragma omp parallel num_threads(ncores)
{
  dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
  lrng.jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps 
  int iter = 0;
#pragma omp for
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      iter = i * n + j;
      // p can be a function of i and j
      boost::random::binomial_distribution<int> dist_binomial(1,p);
      auto gen_bernoulli = std::bind(dist_binomial, std::ref(lrng));
      boost::random::normal_distribution<int> dist_normal1(2.0,1.0);
      auto gen_normal1 = std::bind(dist_normal1, std::ref(lrng));
      boost::random::normal_distribution<int> dist_normal2(4.0,3.0);
      auto gen_normal2 = std::bind(dist_normal2, std::ref(lrng));
      out(iter,0) = gen_bernoulli();
      out(iter,1) = gen_normal1();
      out(iter,2) = gen_normal2();
    }
  }
}
// ok to use rng here
return out;
}

/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/

When I try to run it Rstudio crashes. However, when I change the code like follows it does work:

// [[Rcpp::depends(dqrng, BH, RcppArmadillo)]]
#include <RcppArmadillo.h>
#include <boost/random/binomial_distribution.hpp>
#include <xoshiro.h>
#include <dqrng_distribution.h>
// [[Rcpp::plugins(openmp)]]
#include <omp.h>

// [[Rcpp::plugins(cpp11)]]

// [[Rcpp::export]]
arma::mat parallel_random_matrix(int n, int m, int ncores, double p=0.5) {
  dqrng::xoshiro256plus rng(42);
  arma::mat out(n*m,3);
  // ok to use rng here
#pragma omp parallel num_threads(ncores)
{
  dqrng::xoshiro256plus lrng(rng);      // make thread local copy of rng 
  lrng.jump(omp_get_thread_num() + 1);  // advance rng by 1 ... ncores jumps 
  int iter = 0;
#pragma omp for
  for (int i = 0; i < m; ++i) {
    for (int j = 0; j < n; ++j) {
      iter = i * n + j;
      // p can be a function of i and j
      boost::random::binomial_distribution<int> dist_binomial(1,p);
      auto gen_bernoulli = std::bind(dist_binomial, std::ref(lrng));
      boost::random::normal_distribution<int> dist_normal1(2.0,1.0);
      auto gen_normal1 = std::bind(dist_normal1, std::ref(lrng));
      boost::random::normal_distribution<int> dist_normal2(4.0,3.0);
      auto gen_normal2 = std::bind(dist_normal2, std::ref(lrng));
      out(iter,0) = gen_bernoulli();
      out(iter,1) = 2.0;//gen_normal1();
      out(iter,2) = 3.0;//gen_normal2();
    }
  }
}
// ok to use rng here
return out;
}

/*** R
parallel_random_matrix(5, 5, 4, 0.75)
*/

What am I doing wrong?

Upvotes: 0

Views: 275

Answers (1)

Ralf Stubner
Ralf Stubner

Reputation: 26823

Here lies the problem:

  boost::random::normal_distribution<int> dist_normal1(2.0,1.0);
                                     ^^^

This distribution is meant for real types, not integral types, c.f. https://www.boost.org/doc/libs/1_69_0/doc/html/boost/random/normal_distribution.html. Correct would be

  boost::random::normal_distribution<double> dist_normal1(2.0,1.0);

Upvotes: 2

Related Questions