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. <mat...@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@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@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
