[Rcpp-devel] Best way to compute M'M if M is triangular using RcppArmadillo

Conrad S conradsand.arma at gmail.com
Wed Apr 17 16:06:49 CEST 2013


On Wed, Apr 17, 2013 at 10:24 PM, F.Tusell <fernando.tusell at ehu.es> wrote:
> SEXP Prod1(SEXP M_) {
>   const mat M = as<mat>(M_) ;
>   mat prodM = M ;
>   prodM = trans(M) * M  ;
>   return(wrap(prodM)) ;
> }

Why are you making a copy of M (as prodM) before doing the multiplication ?
Why not just directly do "mat prodM = M.t() * M" ?


> SEXP Prod2(SEXP M_) {
>   const mat M = as<mat>(M_) ;
>   mat prodM = M ;
>   prodM = trimatu(M).t() * trimatu(M)  ;
>   return(wrap(prodM)) ;
> }
> ...
> I expected that telling Armadillo that M is upper triangular would
> make for a fastest product, and the contrary seems true;

trimatu() and trimatl() are currently only used by inv() and solve()
in Armadillo.  Other expressions (such as multiplication) haven't been
optimised yet to take advantage of the extra information on matrix
structure; the trimatu()/trimatl() functions end up creating a
triangular matrix before multiplication (ie. partial copy).


More information about the Rcpp-devel mailing list