2

The following portion of the function rl_compute compiles fine but causes R to crash when I run it. Not sure why. It works with 1 thread.

// [[Rcpp::export]]
void rl_compute(NumericVector covparms, 
        NumericMatrix locs, 
        NumericMatrix NNarray,
        NumericVector& y, 
        NumericMatrix X)
{
    
    int n = y.length();
    int m = NNarray.ncol();
    int p = X.ncol();
    int nparms = covparms.length();
    int dim = locs.ncol();

#pragma omp parallel num_threads(2)\
    {
    int id = omp_get_thread_num();
    for(int i=id; i<n; i= i+2 ){

    int bsize = std::min(i+1,m);

    NumericMatrix locsub(bsize, dim);
    arma::mat X0( bsize, p );
    arma::vec ysub;
    NumericVector idx = na_omit(NNarray.row(i));
    ...}
    }

}
yf297
  • 21
  • 2
  • 5
    Welcome to StackOverflow. For the particular question, read for example the vignettes and examples for RcppParallel that say in no uncertain terms that you can *never ever* have R objects within the parallel sections. So no `NumericMatrix`, `NumericVector` and no Armadillo objects created from R (via Rcpp) via the default (advanced) constructor using zero-copy. – Dirk Eddelbuettel Jul 12 '20 at 21:20
  • @DirkEddelbuettel got it thanks – yf297 Jul 12 '20 at 21:34
  • 3
    So here `X0` and `ysub` should be fine. As could `locsub` be as an `arma::mat`. But you need to create `idx` "freshly" with a distinct copy (so Armadillo docs, or just copy explicitly). Then try again. That may do the trick. If not, well, parallel code debugging is always twice as much fun :-/ – Dirk Eddelbuettel Jul 12 '20 at 21:44
  • A side question: Any particular reason for that fact that you don't use `omp parallel for schedule(static,1) num_threads(2)` but instead reimplement it using `for (int i = omp_get_thread_num(); i < n; i += 2)`? – Hristo Iliev Jul 13 '20 at 10:54

0 Answers0