I'm using Rcpp with OpenMP for parallel computation in C++. I read Is accessing/setting elements in a Rcpp::List thread safe? and it seems that using Rcpp inside a threaded code block is dangerous. But it is quite common to read or write on Rcpp NumericVectors in a for loop such as https://gallery.rcpp.org/articles/hierarchical-risk-parity/.
Below is what I wrote. The task is straightforward. I first declare a NumericVector with a large size n outside a (parallel) for loop. I do some calculation and fill each element of the NumericVector using index i. Once done, I compute a desired quantile of it.
#include <Rcpp.h>
#include <omp.h>
// [[Rcpp::export]]
double foo(const int n, const int nthreads) {
Rcpp::NumericVector out(n);
#pragma omp parallel for num_threads(nthreads)
for (int i = 0; i < n; ++i) {
out[i] = doSomething();
}
Rcpp::Environment stats("package:stats");
Rcpp::Function quantile = stats["quantile"];
return Rcpp::as<double>(quantile(out, Rcpp::Named("probs") = 0.95));
}
I have two questions.
Is there any reason to use other containers (e.g.
std::vector<double>) instead ofNumericVectorfor thread safety? I tested the above code withstd::vector<double>and got identical result. I heard that conversion between the two involves deep copy. If I usestd::vector<double>, then does deep copying occur whenquntile()is called?If using
NumericVectoris safe, would it be better to just returnNumericVectorand applyquntile()in R in terms of performance?