[Rcpp-devel] Parallel random numbers using Rcpp and OpenMP

[Romain Francois]

Hi, 

If you can assume c++11, you might bot need openmp, as you can just use standard support for threads, etc ... 

Although the compiler suite used on windows (if you care about that) does not support c++11 threads. This might be available later. 

Romain


Le 28 avr. 2014 à 12:51, Matteo Fasiolo <matteo.fasiolo at gmail.com> a écrit :

> 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
> 
> _______________________________________________
> 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/9880b658/attachment.html>