This document provides some comparisons of three methods to call C from R: .C, .Call, and Rcpp. Here are descriptions for these methods:

rowSums example

In this example I want to compute the row sums of a double matrix, mimicing the rowSums function in R.

Here is the C function myRowSum.c for using .C and .Call. There are two functions: myRowSum_C is to use .C, and myRowSum_Call uses .Call.

#include <stdio.h>
#include <R.h>
#include <Rinternals.h>

void myRowSum_C(double *x, int *nrowx, int *ncolx, double *result) {
        int i, j;

        for(i=0; i<nrowx[0]; i++)
                for(j=0; j<ncolx[0]; j++)
                        result[i] = result[i] + x[j*nrowx[0]+i];
}

SEXP myRowSum_Call(SEXP X) {
        int i, j, nrowx, ncolx;
        double *x_ptr;

        /* extract the "dim" attribute */
        SEXP dims;
        dims = getAttrib( X, R_DimSymbol ) ;
        nrowx = INTEGER(dims)[0];
        ncolx = INTEGER(dims)[1];
        x_ptr = REAL(X);

        /* allocate memory for the result vector */
        SEXP result = PROTECT( allocVector(REALSXP, nrowx) );
        double *result_ptr = REAL(result);

        for(i=0; i<nrowx; i++)  {
                result_ptr[i] = 0.0;
                for(j=0; j<ncolx; j++)
                        result_ptr[i] = result_ptr[i] + x_ptr[j*nrowx+i];
        }
        UNPROTECT(1);
        return(result);
}

I compiled the above file outside R using R CMD SHLIB myRowSum.c, which generates share library myRowSum.so.

Here is the cpp function myRowSum_Rcpp.cpp for using Rcpp. This function will be compiled in R using sourceCpp R funtion.

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]                                                                                 
NumericVector myRowSum_Rcpp( NumericMatrix x ) {
    int nrow = x.nrow(), ncol = x.ncol();
    int i, j;
    double tmp;
    NumericVector result(nrow);

    for(i=0; i<nrow; i++) {
        tmp = 0.0;
        for(j=0; j<ncol; j++)
            tmp = tmp + x(i,j);
        result[i] = tmp;
    }
    return(result);
}

Below are the R codes for comparison. 1. First load in necessary packages, shared libray, and compile the Rcpp code.

library(Rcpp)
library(microbenchmark)
dyn.load("myRowSum.so")
sourceCpp("myRowSum_Rcpp.cpp")
  1. I write two R functions to wrap myRowSum_C and myRowSum_Call functions from C.
myRowSum_C <- function(X) {
    .C("myRowSum_C", X, as.integer(nrow(X)), as.integer(ncol(X)), rep(0,nrow(X)))[[4]]
}

myRowSum_Call <- function(X) {
    .Call("myRowSum_Call", X)
}
  1. I create a random matrix with 1 million rows and 10 columns, and benchmark the three functions. As a baseline, the rowSums R function is also being compared.
X = matrix(rnorm(1e7,sd=5), ncol=10)
microbenchmark(myRowSum_C(X), myRowSum_Call(X), myRowSum_Rcpp(X), rowSums(X), times=100)
## Unit: milliseconds
##              expr       min        lq      mean    median        uq
##     myRowSum_C(X) 34.848010 40.890005 50.549986 48.308802 54.529563
##  myRowSum_Call(X)  8.314495  9.923796 12.303820 10.662178 11.869721
##  myRowSum_Rcpp(X)  6.197581  7.396726  9.964185  8.096877  9.242808
##        rowSums(X) 21.238879 24.793712 27.591730 26.218080 28.845653
##       max neval cld
##  98.69287   100   c
##  40.09937   100 a  
##  66.39647   100 a  
##  52.33727   100  b

It can be seen that on average, myRowSum_Call takes about 12 milliseconds, myRowSum_Rcpp takes around 10 milliseconds, and myRowSum_C is in the slowest taking 50 milliseconds. As a benchmark, the rowSums function in R takes ~29 milliseconds. Of course that function has some other error checking machenism, which slows it down. These numbers will vary from run to run, but it appears Rcpp is the fastest, followed by .Call. .C is the slowest by a large margin. The discrepency is so big that it is suprising.

A linear kernel smoothing example

The above example only tests the opertaions on loop and sum. Now I want to throw in multiplication. In this example, I implement the linear kernel smoothing procedure, mimicing the filter function in R. This time I will just compare .Call and Rcpp.

Here is the C funnction using .Call myFilter_Call.c:

#include <R.h>
#include <Rinternals.h>

SEXP myFilter_Call( SEXP x, SEXP f ) {
    int nx = length(x), nf = length(f);
    int i, j;
    double tmp;

    SEXP result = PROTECT( allocVector(REALSXP, nx) );
    double *result_ptr = REAL(result);
    double *x_ptr = REAL(x), *f_ptr = REAL(f);

    int half_nf = (nf+1)/2;
    for(i=half_nf-1; i<nx - half_nf; i++) {
        tmp = 0.0;
        for(j=0; j<nf; j++) {
            tmp += x_ptr[i-half_nf+1+j] * f_ptr[j];
        }
        result_ptr[i] = tmp;
    }
    UNPROTECT(1);
    return(result);
}

Here is the cpp funnction using Rcpp myFilter_Rcpp.cpp:

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector myFilter_Rcpp( NumericVector x, NumericVector f ) {
    int nx = x.length(), nf = f.length();
    int i, j;
    double tmp;
    NumericVector result(nx);

    int half_nf = (nf+1)/2;
    for(i=half_nf-1; i<nx - half_nf; i++) {
        tmp = 0.0;
        for(j=0; j<nf; j++) {
            tmp += x[i-half_nf+1+j]*f[j];
        }
        result[i] = tmp;
    }
    return(result);
}

I then use the following R code to do the comparison. Here I use a Gaussian kernel with length 9 to smooth 1 million data points.

library(Rcpp)
library(microbenchmark)
sourceCpp("myFilter_Rcpp.cpp")

dyn.load("myFilter_Call.so")
myFilter_Call <- function(X, f) {
    .Call("myFilter_Call", X, f)
}

X = rnorm(1e6)
f = dnorm(qnorm(seq(0.1, 0.9, by=0.1)))
f = f/sum(f) ## Gaussian kernel       
microbenchmark( myFilter_Call(X,f), myFilter_Rcpp(X,f), times=100)
## Unit: milliseconds
##                 expr      min       lq      mean   median       uq
##  myFilter_Call(X, f) 5.382100 6.465250  9.372279 8.658638  9.99313
##  myFilter_Rcpp(X, f) 6.311626 7.809008 10.110947 9.800286 10.41113
##       max neval cld
##  38.76324   100   a
##  34.25049   100   a

Again, Rcpp and .Call have comparable performances.

Conclusion

Based on the comparisons above, .Call as an old-fashion way to interface with C code provides similar performance as Rcpp. However, Rcpp has the easiest interface, whereas .Call requires a little bit more programming in manipulating the input/output objects. Thus, Rcpp seems to be the best interface to C from R.