On Thu, Jul 14, 2011 at 8:21 AM, Alireza Mahani <alireza.s.mah...@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@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > -- Gabriel Becker Graduate Student Statistics Department University of California, Davis
______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel