[Rd] Performance of .C and .Call functions vs. native R code

[Alireza Mahani]

(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)

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 vectors 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.