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

[Dirk Eddelbuettel]

Matteo, 

On 28 April 2014 at 01:30, Matteo Fasiolo wrote:
| Dear Rcpp developers,
| 
|  I am working on a small package (here) that should provide efficient tools 
| for the multivariate normal distribution using Rcpp/RcppArmadillo and OpenMP.
|  
| Creating functions that evaluate the multivariate normal density efficiently
| was fairly straightforward, 

Yes, and there is also a Rcpp Gallery post on this topic:
   http://gallery.rcpp.org/articles/dmvnorm_arma/

| but I am struggling with parallel random number
| generation with OpenMP.
| 
| 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.

You cannot 'touch R' at all if you are in an OpenMP context as R is nowhere
near thread safe.  Use R as a 'shell' to launch a job, keep data in C++
containers and go crazy in OpenMP. Then collect results and report back to R.
| 
| Is anybody aware of any package providing a parallel C++ rng which my
| package might link to? I have read this post about Tina's rng, which
| seems to work with OpenMP and Rcpp. How hard would it be to have
| such a library included in my package (if at all possible)?

A number of people have worked on this. I used this to try to make the 'port'
of DEoptim to C++ / Rcpp work in parallel.  I got it partially working but
didn't get to a point where I liked the reproducibility.

Around that time, Petr Savicky had already done some work, and I presume
continued with that.  See his page on this research which you can find at
http://www2.cs.cas.cz/~savicky/randomNumbers/ --  I not sure if he has code
at R-Forge or GitHub.  From a glance at the page, it looks like he as a
package for the SFMT generator (by the Mersenne Twister generators) for
OpenMP.

And at the same time Karl Forner started to some work, and also exchanged
notes with Petr.  Karl does have a repo at R-Forge, it has several related
package, and he also had done work on SFMT via OpenMP (and Rcpp). See
eg here: https://r-forge.r-project.org/scm/viewvc.php/pkg/?root=gwas-bin-tests

I am going to CC Petr and Karl. It would be great to settle on a reference
implementation or two, and I am sure between the two of them we will get
further pointers as to what works.

Cheers, Dirk

| 
| Sorry for the possibly silly questions and thanks for any suggestion.
| 
| Matteo
| 
| ----------------------------------------------------------------------
| _______________________________________________
| 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
-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com