(This article was first published on

**Rcpp Gallery**, and kindly contributed to R-bloggers)We want to do matrix multiplication for 3 cases:

- dense times dense
- sparse times dense for sparse matrices of class
`dgCMatrix`

- sparse times dense for sparse matrices of class
`indMatrix`

,

using R’s `Matrix`

package for sparse matrices in R and

`RcppArmadillo`

for C++ linear algebra:

```
// [[Rcpp::depends(RcppArmadillo)]]
#include
```
using namespace Rcpp ;
arma::mat matmult_sp(const arma::sp_mat X, const arma::mat Y){
arma::mat ret = X * Y;
return ret;
};
arma::mat matmult_dense(const arma::mat X, const arma::mat Y){
arma::mat ret = X * Y;
return ret;
};
arma::mat matmult_ind(const SEXP Xr, const arma::mat Y){
// pre-multiplication with index matrix is a permutation of Y's rows:
arma::uvec perm = as<S4>(Xr).slot("perm");
arma::mat ret = Y.rows(perm - 1);
return ret;
};
//[[Rcpp::export]]
arma::mat matmult_cpp(SEXP Xr, const arma::mat Y) {
if (Rf_isS4(Xr)) {
if(Rf_inherits(Xr, "dgCMatrix")) {
return matmult_sp(as<arma::sp_mat>(Xr), Y) ;
} ;
if(Rf_inherits(Xr, "indMatrix")) {
return matmult_ind(Xr, Y) ;
} ;
stop("unknown class of Xr") ;
} else {
return matmult_dense(as<arma::mat>(Xr), Y) ;
}
}

**Set up test cases:**

`library(Matrix)`

Loading required package: methods

```
library(rbenchmark)
set.seed(12211212)
n <- 1000
d <- 50
p <- 30
X <- matrix(rnorm(n*d), n, d)
X_sp <- as(diag(n)[,1:d], "dgCMatrix")
X_ind <- as(sample(1:d, n, rep=TRUE), "indMatrix")
Y <- matrix(1:(d*p), d, p)
```

**Check exception handling:**

`matmult_cpp(as(X_ind, "ngTMatrix"), Y)`

Error: unknown class of Xr

**Dense times dense:**

`all.equal(X%*%Y, matmult_cpp(X, Y))`

[1] TRUE

```
benchmark(X%*%Y,
matmult_cpp(X, Y),
replications=100)[,1:4]
```

test replications elapsed relative 2 matmult_cpp(X, Y) 100 0.086 1.00 1 X %*% Y 100 0.098 1.14

`dgCMatrix`

times dense:

```
all.equal(as(X_sp%*%Y, "matrix"), matmult_cpp(X_sp, Y),
check.attributes = FALSE)
```

[1] TRUE

```
benchmark(X_sp%*%Y,
matmult_cpp(X_sp, Y),
replications=100)[,1:4]
```

test replications elapsed relative 2 matmult_cpp(X_sp, Y) 100 0.009 1.000 1 X_sp %*% Y 100 0.025 2.778

`indMatrix`

times dense:

`all.equal(X_ind%*%Y, matmult_cpp(X_ind, Y))`

[1] TRUE

```
benchmark(X_ind%*%Y,
matmult_cpp(X_ind, Y),
replications=100)[,1:4]
```

test replications elapsed relative 2 matmult_cpp(X_ind, Y) 100 0.013 1.000 1 X_ind %*% Y 100 0.025 1.923

Based on this Q&A on StackOverflow,

thanks again to Kevin Ushey for his helpful comment.

To

**leave a comment**for the author, please follow the link and comment on their blog:**Rcpp Gallery**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...