From: Gabriel Becker <gmbecker_at_ucdavis.edu>

Date: Thu, 14 Jul 2011 08:59:49 -0700

Date: Thu, 14 Jul 2011 08:59:49 -0700

On Thu, Jul 14, 2011 at 8:21 AM, Alireza Mahani <alireza.s.mahani_at_gmail.com>wrote:

> (I am using a LINUX machine)

*>
**> Jeff,
**>
**> In creating reproducible results, I 'partially' answered my question. I
**> have
**> attached two scripts, 'mvMultiply.r' and 'mvMultiply.cc'. Please copy both
**> files into your chosen directory, then run 'Rscript mvMultiply.r' in that
**> directory while changing the two boolean parameters 'INCLUDE_DATAPREP' and
**> 'ROWMAJOR' to all four permutations. (The variable 'diffVec' is there to
**> verify that the two methods produce the same output vector.)
**>
**> Below are the results that I get, along with discussion (tR and tCall are
**> in
**> sec):
**>
**> INCLUDE_DATAPREP,ROWMAJOR,tR,tCall
**> F,F,13.536,13.875
**> F,T,13.824,14.299
**> T,F,13.688,18.167
**> T,T,13.982,30.730
**>
**> Interpretation: The execution time for the .Call line is nearly identical
**> to
**> the call to R operator '%*%'. Two data preparation lines for the .Call
**> method create the overhead:
**>
**> A <- t(A) (~12sec, or 12usec per call)
**> dim(A) <- dim(A)[1] * dim(A)[2] (~4sec, or 4usec per call)
**>
**>
*

AFAIK R stores matrices as vectors internally anyway and the dims just tell it the position of the various elements, so I'm not sure that second line is needed at all.

I have attached a tiny piece of c code which verifies this.

The output I get from that is:

*> dyn.load("/home/gmbecker/gabe/matvectest.so")
**> vec = 1.1:8.1
**> mat = matrix(vec, ncol = 4)
**> .Call("R_MatVecTest", vec, mat, 8L)
*

**[1] TRUE
**
Note if you create the matrix with byrow=TRUE they may not be the same.

Hope that helps,

Gabe

> While the first line can be avoided by providing options in c++ function

*> (as
**> is done in the BLAS API), I wonder if the second line can be avoided, aside
**> from the obvious option of rewriting the R scripts to use evectors instead
**> of
**> matrices. But this defies one of the primary advantages of using R, which
**> is
**> succinct, high-level coding. In particular, if one has several matrices as
**> input into a .Call function, then the overhead from matrix-to-vector
**> transformations can add up. To summarize, my questions are:
**>
**> 1- Do the above results seem reasonable to you? Is there a similar penalty
**> in R's '%*%' operator for transforming matrices to vectors before calling
**> BLAS functions?
**> 2- Are there techniques for reducing the above overhead for developers
**> looking to augment their R code with compiled code?
**>
**> Regards,
**> Alireza
**>
**> ---------------------------------------
**> # mvMultiply.r
**> # comparing performance of matrix multiplication in R (using '%*%'
**> operator)
**> vs. calling compiled code (using .Call function)
**> # y [m x 1] = A [m x n] %*% x [n x 1]
**>
**> rm(list = ls())
**> system("R CMD SHLIB mvMultiply.cc")
**> dyn.load("mvMultiply.so")
**>
**> INCLUDE_DATAPREP <- F
**> ROWMAJOR <- F #indicates whether the c++ function treats A in a row-major
**> or
**> column-major fashion
**>
**> m <- 100
**> n <- 10
**> N <- 1000000
**>
**> diffVec <- array(0, dim = N)
**> tR <- 0.0
**> tCall <- 0.0
**> for (i in 1:N) {
**>
**> A <- runif(m*n); dim(A) <- c(m,n)
**> x <- runif(n)
**>
**> t1 <- proc.time()[3]
**> y1 <- A %*% x
**> tR <- tR + proc.time()[3] - t1
**>
**> if (INCLUDE_DATAPREP) {
**> t1 <- proc.time()[3]
**> }
**> if (ROWMAJOR) { #since R will convert matrix to vector in a
**> column-major
**> fashion, if the c++ function expects a row-major format, we need to
**> transpose A before converting it to a vector
**> A <- t(A)
**> }
**> dim(A) <- dim(A)[1] * dim(A)[2]
**> if (!INCLUDE_DATAPREP) {
**> t1 <- proc.time()[3]
**> }
**> y2 <- .Call("matvecMultiply", as.double(A), as.double(x),
**> as.logical(c(ROWMAJOR)))
**> tCall <- tCall + proc.time()[3] - t1
**>
**> diffVec[i] <- max(abs(y2 - y1))
**> }
**> cat("Data prep time for '.Call' included: ", INCLUDE_DATAPREP, "\n")
**> cat("C++ function expects row-major matrix: ", ROWMAJOR, "\n")
**> cat("Time - Using '%*%' operator in R: ", tR, "sec\n")
**> cat("Time - Using '.Call' function: ", tCall, "sec\n")
**> cat("Maximum difference between methods: ", max(diffVec), "\n")
**>
**> dyn.unload("mvMultiply.so")
**> ---------------------------------------
**> # mvMultiply.cc
**> #include <Rinternals.h>
**> #include <R.h>
**>
**> extern "C" {
**>
**> SEXP matvecMultiply(SEXP A, SEXP x, SEXP rowmajor) {
**> double *rA = REAL(A), *rx = REAL(x), *ry;
**> int *rrm = LOGICAL(rowmajor);
**> int n = length(x);
**> int m = length(A) / n;
**> SEXP y;
**> PROTECT(y = allocVector(REALSXP, m));
**> ry = REAL(y);
**> for (int i = 0; i < m; i++) {
**> ry[i] = 0.0;
**> for (int j = 0; j < n; j++) {
**> if (rrm[0] == 1) {
**> ry[i] += rA[i * n + j] * rx[j];
**> } else {
**> ry[i] += rA[j * m + i] * rx[j];
**> }
**> }
**> }
**> UNPROTECT(1);
**> return(y);
**> }
**>
**> }
**>
**>
**> --
**> View this message in context:
**> http://r.789695.n4.nabble.com/Performance-of-C-and-Call-functions-vs-native-R-code-tp3665017p3667896.html
**> Sent from the R devel mailing list archive at Nabble.com.
**>
**> ______________________________________________
**> R-devel_at_r-project.org mailing list
**> https://stat.ethz.ch/mailman/listinfo/r-devel
**>
*

-- Gabriel Becker Graduate Student Statistics Department University of California, DavisReceived on Thu 14 Jul 2011 - 16:09:38 GMT______________________________________________ R-devel_at_r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

*
Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.
Archive generated by hypermail 2.2.0, at Thu 14 Jul 2011 - 16:40:09 GMT.
*

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*