[Rcpp-devel] Best way to create random normal matrix

[Dirk Eddelbuettel]

Julian,

Your prorgam was actually fine as such, but you "used" it wrong -- what comes
back from cxxfunction() is already an R function, see below.

If you use the R RNGs (which is easy from Rcpp) your results are reproducible
with R.  That is all documented so I won't repeat it here.

Dirk

R> library(inline)
R> 
R> src <- '
+     Rcpp::NumericMatrix Xr(Xs);
+     int q = Rcpp::as<int>(ys);
+     int n = Xr.nrow(), k = Xr.ncol();
+     arma::mat X(Xr.begin(), n, k, false);
+     arma::mat G, Y, B;
+     G = arma::randn(n,q);
+     Y = X*G;
+     arma::mat Q, R;
+     arma::qr(Q,R,Y);
+     Q = Q(arma::span::all,arma::span(0,q-1));
+     B = arma::trans(Q)*X;
+     arma::vec d;
+     arma::mat u, v;
+     arma::svd(u, d, v, B );
+     arma::mat U = Q*u;
+     return Rcpp::List::create(Rcpp::Named("d")=d,
+                               Rcpp::Named("U")=U,
+                               Rcpp::Named("V")=v);'
R> 
R> rsvd <- cxxfunction(signature(Xs="numeric", ys="integer"), body=src, plugin= "RcppArmadillo")
R> rsvd(matrix(1:4,2,2), 2)
$d
         [,1]
[1,] 5.464986
[2,] 0.365966

$U
         [,1]      [,2]
[1,] 0.576048 -0.817416
[2,] 0.817416  0.576048

$V
         [,1]      [,2]
[1,] 0.404554  0.914514
[2,] 0.914514 -0.404554

R> 


-- 
Dirk Eddelbuettel | edd at debian.org | http://dirk.eddelbuettel.com