[Rcpp-devel] Parallel random numbers using Rcpp and OpenMP
Hi all,
many thanks for your detailed replies, I see that there are many options
out there.
As a first step I think I will try out the C++11 option that Matt
suggested, as it seems
relatively simple. In particular I might do something like:
#pragma omp parallel
{
std::mt19937_64 engine(
static_cast<uint64_t>(seeds[omp_get_thread_num()] ) );
std::uniform_real_distribution<double> zeroToOne(0.0, 1.0);
#pragma omp for
for (int i = 0; i < 1000; i++) {
double a = zeroToOne(engine);
}
}
were seeds is a NumericVector coming from R. That's probably not ideal, but
it might
give reasonable and reproducible results.
Thanks,
Matteo
On Mon, Apr 28, 2014 at 9:39 AM, Matt D. <matdzb at gmail.com> wrote:
On 4/28/2014 01:30, Matteo Fasiolo wrote:
[...] As I understand, doing things such as: NumericMatrix out(n, d); #pragma omp for schedule(static) for(int kk = 0; kk < d; kk++) out( _, kk) = rnorm(n); is not going to work, because rnorm() is not thread safe (in fact this code crashed R). On the other hand R level parallelism using clusterApply, mclapply etc appears to be too slow to be of any use for this purpose. Is anybody aware of any package providing a parallel C++ rng which my package might link to? Your first choice can be to go with the C++ Standard Library -- header
<random> -- (importantly) using one object per thread: http://stackoverflow.com/questions/8813592/c11-thread- safety-of-random-number-generators See: http://en.cppreference.com/w/cpp/numeric/random http://isocpp.org/blog/2013/03/n3551-random-number-generation Since you're using OpenMP, you can get the distinct thread ID by using the `omp_get_thread_num` function: http://stackoverflow.com/questions/15918758/how-to- make-each-thread-use-its-own-rng-in-c11 Note that if you're doing large scale parallel statistics (say, Monte Carlo) you may also want a PRNG with statistical properties which don't deteriorate (e.g., no spurious correlation, let alone overlapping) for very large samples; preferably, also one that doesn't maintain any kind of global state (so-called "stateless") is going to be easier to use, too (nothing to synchronize/lock across threads). A relatively recent work that's particularly relevant is Random123: https://www.deshawresearch.com/resources_random123.html http://www.thesalmons.org/john/random123/ MIT-licensed C++ version is available here: http://www.sitmo.com/article/parallel-random-number-generator-in-c/ With a simple example: http://www.sitmo.com/article/multi-threaded-random-number- generation-in-c11/ // The author is currently working on getting it into Boost.Random; discussion: http://www.wilmott.com/messageview.cfm?catid=44& threadid=95882&STARTPAGE=2#710955 HTH, Best, Matt
_______________________________________________ Rcpp-devel mailing list Rcpp-devel at lists.r-forge.r-project.org https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-------------- next part -------------- An HTML attachment was scrubbed... URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20140428/ea6298f3/attachment.html>