[R-pkg-devel] Cannot create C code with acceptable performance with respect to internal R command.

[Dirk Eddelbuettel]

Luc,

As Tomas mentioned, matrix-multiplication can take advantage of multiple
threads, and the 'text book' nexted loops do not do that.  Now, one
alternative that appeals a lot to me is to farm out to Armadillo which also
calls LAPACK for you (as R does). And via RcppArmadillo, the setup becomes a
one-liner with the expression 'mat1 * mat2' where '*' is overloaded
appropriately (as is matrix multiplication '%*%' in R).  I include your
example as self-contained and reproducible script below, on my not-so-recent
machine with twelve cores I get

$ Rscript luc.r 
Unit: microseconds
 expr       min        lq     mean    median       uq      max neval cld
    C 29010.538 39242.004 47948.98 50930.500 52715.30 81668.53   100  a 
    R   685.658   800.653  1984.17  1129.754  2719.88  8420.66   100   b
  Cpp   401.182   444.164  1775.03   651.023  1656.24 30369.15   100   b
$ 

but what really shines (in my eyes) is that a function

    arma::mat cppprod(const arma::mat& m1, const arma::mat& m2) {
        return m1 * m2;
    }

gets set-up for you with no worries whatsoever and outscores the R
version. (And if you look into the Rcpp docs you can learn to make this a
little faster still but skipping a (generally recommended !!) handshake with
RNG status etc).

But different strokes for different folks, not everybody likes C++ (which is
both perfectly find and also includes Tomas who saw fit to rail against it
yesterday regarding its compile times which can both tweaked and are also
worse still in some other popular languages) but I digress ...

Hope this helps, Dirk


ccode <- r"(
SEXP u1 = Rf_getAttrib(mat1, R_DimSymbol);
int m1 = INTEGER(u1)[0];
int n1 = INTEGER(u1)[1];
SEXP u2 = Rf_getAttrib(mat2, R_DimSymbol);
int m2 = INTEGER(u2)[0];
int n2 = INTEGER(u2)[1];
if (n1 != m2) Rf_error("matrices not conforming");
SEXP retval = PROTECT(Rf_allocMatrix(REALSXP, m1, n2));
double* left = REAL(mat1);
double* right = REAL(mat2);
double* ret = REAL(retval);
double werk = 0.0;
for (int j = 0; j < n2; j++) {
  for (int i = 0; i < m1; i++) {
     werk = 0.0;
     for (int k = 0; k < n1; k++)
       werk += (left[i + m1 * k] * right[k + m2 * j]);
     ret[j * m1 + i] = werk;
  }
}
UNPROTECT(1);
return retval;
)"
cprod <- inline::cfunction(sig=signature(mat1="numeric", mat2="numeric"), body=ccode, language="C")

Rcpp::cppFunction("arma::mat cppprod(const arma::mat& m1, const arma::mat& m2) { return m1 * m2; }", depends="RcppArmadillo")

set.seed(123)
m1 <- matrix(rnorm(300000), nrow = 60)
m2 <- matrix(rnorm(300000), ncol = 60)
print(microbenchmark::microbenchmark(C = cprod(m1, m2),
                                     R = m1 %*% m2,
                                     Cpp = cppprod(m1, m2),
                                     times = 100))

-- 
dirk.eddelbuettel.com | @eddelbuettel | edd using debian.org