randu 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_randu(size_t m, size_t n, double min = 0.0, double max = 1.0) {
  static std::random_device rd;
  static std::uniform_real_distribution<double> dist(min, max);

  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_randu_omp(size_t m, size_t n, double min = 0.0, double max = 1.0) {
  blaze::setNumThreads( blaze::getNumThreads() );
  static std::random_device rd;
  static std::uniform_real_distribution<double> dist(min, max);

  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_randu(arma::uword m, arma::uword n, double min = 0.0, double max = 1.0) {
  return arma::randu<arma::mat>(m, n) * (max - min) + min;
}")

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

// [[Rcpp::export]]
Eigen::MatrixXd rcppeigen_randu(std::size_t m, std::size_t n, double min = 0.0, double max = 1.0) {
  Eigen::MatrixXd A = Eigen::MatrixXd::Random(m, n).cwiseAbs() * (max - min);
  A.array() += min;
  return A;
}")

library(microbenchmark)
nn <- 3e3L
microbenchmark(rcppblaze = rcppblaze_randu(nn, nn, -5, 5),
               rcppblaze_omp = rcppblaze_randu_omp(nn, nn, -5, 5),
               rcpparma = rcpparma_randu(nn, nn, -5, 5),
               rcppeigen = rcppeigen_randu(nn, nn, -5, 5),
               r = matrix(runif(nn*nn), nn)* 10 - 5,
               times = 20L)
# Unit: milliseconds
#           expr       min        lq      mean    median        uq       max neval
#      rcppblaze 116.73880 117.48573 124.39618 120.20502 130.55786 144.13320    20
#  rcppblaze_omp  45.54567  49.05225  55.64187  51.51919  64.37526  75.13228    20
#       rcpparma 128.17825 128.81605 138.73946 131.21512 146.11872 200.86437    20
#      rcppeigen 179.56535 183.39844 190.49903 185.19832 197.97525 208.82546    20
#              r 273.87377 285.63416 292.84457 289.33341 297.28060 351.80331    20