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

F.Tusell fernando.tusell at ehu.es
Wed Apr 17 14:24:17 CEST 2013

```Not a question strictly about Rcpp but hope this is a right

I am trying to find out what is the fastest possible way to compute
M'M for M upper triangular. I have coded,

// [[Rcpp::export]]
SEXP Prod1(SEXP M_) {
const mat M = as<mat>(M_) ;
mat prodM = M ;
prodM = trans(M) * M  ;
return(wrap(prodM)) ;
}

// [[Rcpp::export]]
SEXP Prod2(SEXP M_) {
const mat M = as<mat>(M_) ;
mat prodM = M ;
prodM = trimatu(M).t() * trimatu(M)  ;
return(wrap(prodM)) ;
}

// [[Rcpp::export]]
SEXP Prod3(SEXP M_) {
mat M = as<mat>(M_) ;
int d = M.n_rows ;
mat prodM = M ;
double * vM = M.memptr() ;
double * vprodM = prodM.memptr() ;
const double one = 1 ;
const double zero = 0 ;
F77_CALL(dsyrk)("L","T",&d,&d,&one,vM,&d,&zero,vprodM,&d)  ;
return(wrap(prodM)) ;
}

and then tested with:

require(microbenchmark)
sourceCpp("prods.cpp")
d <- 50
a <- chol(crossprod(matrix(rnorm(d*d),d,d)))

m <- microbenchmark(
r1 <- Prod1(a),
r2 <- Prod2(a),
r3 <- Prod3(a),
times=100
)

This is what I get:

> m
Unit: microseconds
expr     min       lq   median       uq      max neval
r1 <- Prod1(a) 138.749 144.2260 148.4815 159.0830 2456.146   100
r2 <- Prod2(a) 296.193 320.0275 329.4770 342.4185 2763.041   100
r3 <- Prod3(a) 132.150 138.2590 140.9270 152.3675  218.719   100

Prod3 using BLAS dsyrk is about as fast as Prod1, using Armadillo.
I expected that telling Armadillo that M is upper triangular would
make for a fastest product, and the contrary seems true; Prod2 takes
about twice as much time as Prod1 and Prod3.

Is there a faster way?

ft.

```