Thanks for the feedback. This should be easy enough to port over to RcppEigen, right?<br><br><div class="gmail_quote">On Tue, Dec 4, 2012 at 6:17 PM, Dirk Eddelbuettel <span dir="ltr"><<a href="mailto:edd@debian.org" target="_blank">edd@debian.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
Julian,<br>
<br>
Your prorgam was actually fine as such, but you "used" it wrong -- what comes<br>
back from cxxfunction() is already an R function, see below.<br>
<br>
If you use the R RNGs (which is easy from Rcpp) your results are reproducible<br>
with R.  That is all documented so I won't repeat it here.<br>
<br>
Dirk<br>
<br>
R> library(inline)<br>
R><br>
R> src <- '<br>
+     Rcpp::NumericMatrix Xr(Xs);<br>
+     int q = Rcpp::as<int>(ys);<br>
+     int n = Xr.nrow(), k = Xr.ncol();<br>
+     arma::mat X(Xr.begin(), n, k, false);<br>
+     arma::mat G, Y, B;<br>
+     G = arma::randn(n,q);<br>
+     Y = X*G;<br>
+     arma::mat Q, R;<br>
+     arma::qr(Q,R,Y);<br>
+     Q = Q(arma::span::all,arma::span(0,q-1));<br>
+     B = arma::trans(Q)*X;<br>
+     arma::vec d;<br>
+     arma::mat u, v;<br>
+     arma::svd(u, d, v, B );<br>
<div class="im">+     arma::mat U = Q*u;<br>
</div>+     return Rcpp::List::create(Rcpp::Named("d")=d,<br>
+                               Rcpp::Named("U")=U,<br>
+                               Rcpp::Named("V")=v);'<br>
R><br>
R> rsvd <- cxxfunction(signature(Xs="numeric", ys="integer"), body=src, plugin= "RcppArmadillo")<br>
R> rsvd(matrix(1:4,2,2), 2)<br>
$d<br>
         [,1]<br>
[1,] 5.464986<br>
[2,] 0.365966<br>
<br>
$U<br>
         [,1]      [,2]<br>
[1,] 0.576048 -0.817416<br>
[2,] 0.817416  0.576048<br>
<br>
$V<br>
         [,1]      [,2]<br>
[1,] 0.404554  0.914514<br>
[2,] 0.914514 -0.404554<br>
<br>
R><br>
<div class="HOEnZb"><div class="h5"><br>
<br>
--<br>
Dirk Eddelbuettel | <a href="mailto:edd@debian.org">edd@debian.org</a> | <a href="http://dirk.eddelbuettel.com" target="_blank">http://dirk.eddelbuettel.com</a><br>
</div></div></blockquote></div><br>