<div dir="ltr">Hi all,<div><br></div><div> many thanks for your detailed replies, I see that there are many options out there.</div><div><br></div><div>As a first step I think I will try out the C++11 option that Matt suggested, as it seems</div>
<div>relatively simple. In particular I might do something like:</div><div><br></div><div><div>#pragma omp parallel</div><div> {</div><div> std::mt19937_64 engine( static_cast<uint64_t>(seeds[omp_get_thread_num()] ) );</div>
<div> std::uniform_real_distribution<double> zeroToOne(0.0, 1.0);</div><div> </div><div> #pragma omp for</div><div> for (int i = 0; i < 1000; i++) {</div><div> double a = zeroToOne(engine);</div>
<div> }</div><div> }</div></div><div><br></div><div>were seeds is a NumericVector coming from R. That's probably not ideal, but it might</div><div>give reasonable and reproducible results.</div><div><br></div>
<div>Thanks,</div><div><br></div><div>Matteo</div><div><br></div></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Mon, Apr 28, 2014 at 9:39 AM, Matt D. <span dir="ltr"><<a href="mailto:matdzb@gmail.com" target="_blank">matdzb@gmail.com</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">On 4/28/2014 01:30, Matteo Fasiolo wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
[...]<div class=""><br>
<br>
As I understand, doing things such as:<br>
<br>
NumericMatrix out(n, d);<br>
#pragma omp for schedule(static)<br>
for(int kk = 0; kk < d; kk++) out( _, kk) = rnorm(n);<br>
<br>
is not going to work, because rnorm() is not thread safe<br>
(in fact this code crashed R). On the other hand R level parallelism<br>
using clusterApply, mclapply etc appears to be too slow to be of any<br>
use for this purpose.<br>
<br>
Is anybody aware of any package providing a parallel C++ rng which my<br>
package might link to?<br>
<br>
</div></blockquote>
Your first choice can be to go with the C++ Standard Library -- header <random> -- (importantly) using one object per thread:<br>
<a href="http://stackoverflow.com/questions/8813592/c11-thread-safety-of-random-number-generators" target="_blank">http://stackoverflow.com/<u></u>questions/8813592/c11-thread-<u></u>safety-of-random-number-<u></u>generators</a><br>
<br>
See:<br>
<a href="http://en.cppreference.com/w/cpp/numeric/random" target="_blank">http://en.cppreference.com/w/<u></u>cpp/numeric/random</a><br>
<a href="http://isocpp.org/blog/2013/03/n3551-random-number-generation" target="_blank">http://isocpp.org/blog/2013/<u></u>03/n3551-random-number-<u></u>generation</a><br>
<br>
Since you're using OpenMP, you can get the distinct thread ID by using the `omp_get_thread_num` function:<br>
<a href="http://stackoverflow.com/questions/15918758/how-to-make-each-thread-use-its-own-rng-in-c11" target="_blank">http://stackoverflow.com/<u></u>questions/15918758/how-to-<u></u>make-each-thread-use-its-own-<u></u>rng-in-c11</a><br>
<br>
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).<br>
<br>
A relatively recent work that's particularly relevant is Random123:<br>
<a href="https://www.deshawresearch.com/resources_random123.html" target="_blank">https://www.deshawresearch.<u></u>com/resources_random123.html</a><br>
<a href="http://www.thesalmons.org/john/random123/" target="_blank">http://www.thesalmons.org/<u></u>john/random123/</a><br>
<br>
MIT-licensed C++ version is available here:<br>
<a href="http://www.sitmo.com/article/parallel-random-number-generator-in-c/" target="_blank">http://www.sitmo.com/article/<u></u>parallel-random-number-<u></u>generator-in-c/</a><br>
<br>
With a simple example:<br>
<a href="http://www.sitmo.com/article/multi-threaded-random-number-generation-in-c11/" target="_blank">http://www.sitmo.com/article/<u></u>multi-threaded-random-number-<u></u>generation-in-c11/</a><br>
<br>
// The author is currently working on getting it into Boost.Random; discussion: <a href="http://www.wilmott.com/messageview.cfm?catid=44&threadid=95882&STARTPAGE=2#710955" target="_blank">http://www.wilmott.com/<u></u>messageview.cfm?catid=44&<u></u>threadid=95882&STARTPAGE=2#<u></u>710955</a><br>
<br>
HTH,<br>
<br>
Best,<br>
<br>
Matt<div class="HOEnZb"><div class="h5"><br>
<br>
______________________________<u></u>_________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org" target="_blank">Rcpp-devel@lists.r-forge.r-<u></u>project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-<u></u>project.org/cgi-bin/mailman/<u></u>listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br></div>