randn matrix in Rcpp

Rcpp::sourceCpp(code = "
// [[Rcpp::depends(RcppBlaze)]]
// [[Rcpp::plugins(cpp11)]]
#include <RcppBlaze.h>
#include <random>
#include <functional>

// [[Rcpp::export]]
blaze::DynamicMatrix<double> rcppblaze_randn(size_t m, size_t n, double mean = 0.0, double stddev = 1.0) {
  static std::random_device rd;
  static std::normal_distribution<double> dist( mean, stddev );

  blaze::DynamicMatrix<double> A(m, n);
  for( size_t i=0UL; i<A.columns(); ++i)
    std::generate(A.begin(i), A.end(i), std::bind(dist, std::default_random_engine( rd() ) ) );
  return A;
}")

Rcpp::sourceCpp(code = "
// [[Rcpp::depends(RcppBlaze)]]
// [[Rcpp::plugins(cpp11,openmp)]]
#include <RcppBlaze.h>
#include <random>
#include <functional>

// [[Rcpp::export]]
blaze::DynamicMatrix<double> rcppblaze_randn_omp(size_t m, size_t n, double mean = 0.0, double stddev = 1.0) {
  blaze::setNumThreads( blaze::getNumThreads() );
  static std::random_device rd;
  static std::normal_distribution<double> dist( mean, stddev );

  blaze::DynamicMatrix<double> A(m, n);
  #pragma omp parallel for
  for( size_t i=0UL; i<A.columns(); ++i)
    std::generate(A.begin(i), A.end(i), std::bind(dist, std::default_random_engine( rd() ) ) );
  return A;
}")

Rcpp::sourceCpp(code = "
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
arma::mat rcpparma_randn(arma::uword m, arma::uword n, double mean = 0.0, double stddev = 1.0) {
  return arma::randn<arma::mat>(m, n) * stddev + mean;
}")

Rcpp::sourceCpp(code = "
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::plugins(cpp11)]]
#include <RcppEigen.h>
#include <random>

inline double randn(double dummy, double mean, double stddev) {
  static std::random_device rd;
  static std::default_random_engine randEng(rd());
  static std::normal_distribution<double> dist(mean, stddev);
  return dist(randEng);
}

// [[Rcpp::export]]
SEXP rcppeigen_randn(std::size_t m, std::size_t n, double mean = 0.0, double stddev = 1.0) {
  auto randn_f = [&mean, &stddev](double x){ return randn(x, mean, stddev); };
  Eigen::MatrixXd A = Eigen::MatrixXd::Zero(m, n).unaryExpr(randn_f);
  return Rcpp::wrap(A);
}")

library(microbenchmark)
nn <- 3e3L
microbenchmark(rcppblaze = rcppblaze_randn(nn, nn, 30, 20),
               rcppblaze_omp = rcppblaze_randn_omp(nn, nn, 30, 20),
               rcpparma = rcpparma_randn(nn, nn, 30, 20),
               rcppeigen = rcppeigen_randn(nn, nn, 30, 20),
               r = matrix(rnorm(nn*nn), nn)*20 + 30,
               times = 20L)
# Unit: milliseconds
#           expr       min        lq     mean   median       uq      max neval
#      rcppblaze 354.75299 361.32248 371.5106 373.1244 378.2249 393.6559    20
#  rcppblaze_omp  84.60684  97.30579 106.6932 104.9424 117.8100 144.1054    20
#       rcpparma 697.18385 702.32619 714.4662 712.7235 722.5486 750.9504    20
#      rcppeigen 368.96569 373.56575 386.1028 381.8487 395.9508 432.4745    20
#              r 587.96085 595.83358 619.2839 609.2916 624.7770 700.9972    20