Dear all,
When forming the matrix-vector-product one can form the inner products between
As rows and x (mvprod1 below) or form a linear combination of As columns
(mvprod2 below). The difference in computing time is striking (at least for me)
for large matrices and it really illustrates the "caching issue" (I assume):
> microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)
Unit: milliseconds
expr min lq median
uq max neval
mvprod1(A, x) 35.783812 36.310956 36.600186 36.991113 38.922889 100
mvprod2(A, x) 4.408892 4.436882 4.475368 4.568668 5.477176 100
A %*% x 8.514091 8.592230 8.734978 8.811951 9.738653
100
Just wanted to share this!
Cheers
Søren
#include <Rcpp.h>
using namespace Rcpp;
//[[Rcpp::export]]
NumericVector mvprod1 (NumericMatrix A, NumericVector x){
int nrA=A.nrow(), ncA=A.ncol(), i, j;
NumericVector y(nrA);
for (i=0; i<nrA; ++i){
for (j=0; j<ncA; ++j){
y[i] += A(i,j)*x[j];
}
}
return y;
}
//[[Rcpp::export]]
NumericVector mvprod2 (NumericMatrix A, NumericVector x){
int nrA=A.nrow(), ncA=A.ncol(), i, j;
NumericVector y(nrA);
for (j=0; j<ncA; ++j){
for (i=0; i<nrA; ++i){
y[i] += A(i,j)*x[j];
}
}
return y;
}
/*** R
d <- 5
A <- matrix(1.0*(1:d^2),nrow=d)
x <- 1.0*(1:d)
A%*%x
mvprod1(A,x)
mvprod2(A,x)
d <- 2000
A <- matrix(1.0*(1:d^2),nrow=d)
x <- 1.0*(1:d)
library(microbenchmark)
microbenchmark(mvprod1(A,x), mvprod2(A,x), A%*%x)
*/
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
