<div dir="ltr"><div><div><div>I still don't get why you are copying the whole object.<br>But if you really want to initialize an array of the desired size, you can use arma::copy_size<br></div>In you case it will be something like<br>
</div>........<br>mat prodM;<br></div>prodM.copy_size(M);<br>.........<br></div><div class="gmail_extra"><br><br><div class="gmail_quote">On Wed, Apr 17, 2013 at 3:24 PM, F.Tusell <span dir="ltr"><<a href="mailto:fernando.tusell@ehu.es" target="_blank">fernando.tusell@ehu.es</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">El mié, 17-04-2013 a las 14:01 +0000, Smith, Dale (Norcross) escribió:<br>
<div class="im">> I agree with Dirk. Take a look at the arma::mat constructor in fastLm. If you are not modifying M, use the raw memory constructor.<br>
><br>
> <a href="http://gallery.rcpp.org/articles/fast-linear-model-with-armadillo/" target="_blank">http://gallery.rcpp.org/articles/fast-linear-model-with-armadillo/</a><br>
><br>
> I don't see any other way to avoid using Rcpp::wrap in the return statement. Note there may be other places you copy memory, for instance, in Prod1 you write<br>
><br>
> |   mat prodM = M ;<br>
> |   prodM = trans(M) * M  ;<br>
><br>
<br>
</div>True. But in Prod3 I needed a pre-existent target array to store the<br>
result, so I defined it also in Prod1 and Prod2 to make things alike.<br>
Taking a copy of M is indeed a very clumsy way to initialize an array<br>
of the right size.<br>
<div><div class="h5"><br>
> I think the first statement is not necessary.<br>
><br>
> Dale Smith, Ph.D.<br>
> Senior Financial Quantitative Analyst<br>
> Risk & Compliance<br>
> Fiserv<br>
> Office: 678-375-5315<br>
> <a href="http://www.fiserv.com" target="_blank">www.fiserv.com</a><br>
><br>
><br>
> -----Original Message-----<br>
> From: <a href="mailto:rcpp-devel-bounces@r-forge.wu-wien.ac.at">rcpp-devel-bounces@r-forge.wu-wien.ac.at</a> [mailto:<a href="mailto:rcpp-devel-bounces@r-forge.wu-wien.ac.at">rcpp-devel-bounces@r-forge.wu-wien.ac.at</a>] On Behalf Of Dirk Eddelbuettel<br>

> Sent: Wednesday, April 17, 2013 9:17 AM<br>
> To: <a href="mailto:fernando.tusell@ehu.es">fernando.tusell@ehu.es</a><br>
> Cc: <a href="mailto:rcpp-devel@r-forge.wu-wien.ac.at">rcpp-devel@r-forge.wu-wien.ac.at</a><br>
> Subject: Re: [Rcpp-devel] Best way to compute M'M if M is triangular using RcppArmadillo<br>
><br>
><br>
> On 17 April 2013 at 14:24, F.Tusell wrote:<br>
> | Not a question strictly about Rcpp but hope this is a right place to<br>
> | ask.<br>
> |<br>
> | I am trying to find out what is the fastest possible way to compute<br>
> | M'M for M upper triangular. I have coded,<br>
> |<br>
> | // [[Rcpp::export]]<br>
> | SEXP Prod1(SEXP M_) {<br>
> |   const mat M = as<mat>(M_) ;<br>
> |   mat prodM = M ;<br>
> |   prodM = trans(M) * M  ;<br>
> |   return(wrap(prodM)) ;<br>
> | }<br>
> |<br>
> | // [[Rcpp::export]]<br>
> | SEXP Prod2(SEXP M_) {<br>
> |   const mat M = as<mat>(M_) ;<br>
> |   mat prodM = M ;<br>
> |   prodM = trimatu(M).t() * trimatu(M)  ;<br>
> |   return(wrap(prodM)) ;<br>
> | }<br>
> |<br>
> |<br>
> | // [[Rcpp::export]]<br>
> | SEXP Prod3(SEXP M_) {<br>
> |   mat M = as<mat>(M_) ;<br>
> |   int d = M.n_rows ;<br>
> |   mat prodM = M ;<br>
> |   double * vM = M.memptr() ;<br>
> |   double * vprodM = prodM.memptr() ;<br>
> |   const double one = 1 ;<br>
> |   const double zero = 0 ;<br>
> |   F77_CALL(dsyrk)("L","T",&d,&d,&one,vM,&d,&zero,vprodM,&d)  ;<br>
> |   return(wrap(prodM)) ;<br>
> | }<br>
> |<br>
> | and then tested with:<br>
> |<br>
> | require(RcppArmadillo)<br>
> | require(microbenchmark)<br>
> | sourceCpp("prods.cpp")<br>
> | d <- 50<br>
> | a <- chol(crossprod(matrix(rnorm(d*d),d,d)))<br>
> |<br>
> | m <- microbenchmark(<br>
> |        r1 <- Prod1(a),<br>
> |        r2 <- Prod2(a),<br>
> |        r3 <- Prod3(a),<br>
> |        times=100<br>
> |            )<br>
> |<br>
> | This is what I get:<br>
> |<br>
> | > m<br>
> | Unit: microseconds<br>
> |            expr     min       lq   median       uq      max neval<br>
> |  r1 <- Prod1(a) 138.749 144.2260 148.4815 159.0830 2456.146   100<br>
> |  r2 <- Prod2(a) 296.193 320.0275 329.4770 342.4185 2763.041   100<br>
> |  r3 <- Prod3(a) 132.150 138.2590 140.9270 152.3675  218.719   100<br>
> |<br>
> | Prod3 using BLAS dsyrk is about as fast as Prod1, using Armadillo.<br>
> | I expected that telling Armadillo that M is upper triangular would<br>
> | make for a fastest product, and the contrary seems true; Prod2 takes<br>
> | about twice as much time as Prod1 and Prod3.<br>
> |<br>
> | Is there a faster way?<br>
><br>
> You are going about this the right way by starting with empirics.<br>
><br>
> Now, d=50 is not big so your (relative) savings will be rather small.  Plus, the way you instantiate the Arma object DOES create extra copies on the arma object instantiation [ see discussion this list just last week ] as well as on exit -- which may dominate the timing. Try d=100, 200, ... 500 for comparison.<br>

><br>
> Lastly, you could write prod3 without going to arma.<br>
><br>
> Dirk<br>
<br>
</div></div>I did try larger d, and to my surprise the relative timings stay about<br>
the same. Conrad S provides one piece of information I was missing, that<br>
trimatu() and trimatl() only have (as yet) effect with inv() and<br>
solve(), hence their use in Prod3 is just a waste of time. This explains<br>
the results.<br>
<br>
For the time being, I will stick with trans(M) * M. A function like R's<br>
crossprod() in Armadillo would be great. Thank you all for your help.<br>
<span class="HOEnZb"><font color="#888888"><br>
ft.<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
_______________________________________________<br>
Rcpp-devel mailing list<br>
<a href="mailto:Rcpp-devel@lists.r-forge.r-project.org">Rcpp-devel@lists.r-forge.r-project.org</a><br>
<a href="https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel" target="_blank">https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel</a><br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br><div dir="ltr">Ahmadou H. DICKO<br>statistician and applied economist<br>PhD student in Climate change economics<br>Faculty of economics and managment - Cheikh Anta Diop University<br>
West African Science Service Center on Climate Change and Adaptated Land Use (WASCAL)<br>Center for Development Research (ZEF) - University of Bonn <br>twitter : @dickoah<br>github : <a href="https://github.com/dickoa" target="_blank">github/dickoa</a><br>
tel : +221 33 827 55 16<br>portable: +221 77 123 81 69<br></div>
</div>