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) 33.618934 39.547585 50.615149 46.838690 55.334390
##  myRowSum_Call(X)  8.022451  9.591825 11.692959 10.305188 12.143728
##  myRowSum_Rcpp(X)  5.884247  7.369254  9.580274  8.200964  9.453613
##        rowSums(X) 23.288751 25.007887 29.334679 27.468234 30.897620
##        max neval cld
##  116.05325   100   c
##   31.60609   100 a  
##   25.33447   100 a  
##   58.41985   100  b

It can be seen that on average, myRowSum_Call takes about 13 milliseconds, but myRowSum_Rcpp takes around 10 milliseconds, and myRowSum_C is in the slowest taking 55 milliseconds. As a benchmark, the rowSums function in R takes ~28 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.328045 6.786385  9.661941 9.367134 10.50744
##  myFilter_Rcpp(X, f) 6.422510 7.762425 10.420926 9.694676 10.73079
##       max neval cld
##  43.85907   100   a
##  34.45560   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.