Combining a C Routine with Rcpp

Calling a hand-written C routine from an Rcpp function

programming
cpp
Author

James Balamuta

Published

January 1, 2019

This post is a version of the rcpp-and-c README from the R&D Rcpp collection of compiled-code examples. The repository holds the full, runnable source.

This small example R package demonstrates how to call a hand-written C routine, one that talks to R through R’s C API and the SEXP type, from within a C++ function exposed to R with Rcpp. It is one of the “R&D Rcpp” prototype packages, a collection of minimal examples for integrating compiled code with R. Use it as a template whenever you need to wrap existing C code behind an Rcpp interface.

The cRcpp R package provides an example of implementing and using a C routine with C++ code through Rcpp.

Usage

To install the package, you must first have a compiler on your system that is compatible with R. For help on obtaining a compiler consult either macOS or Windows guides.

With a compiler in hand, one can then install the package from GitHub by:

# install.packages("remotes")
remotes::install_github("coatless-rd-rcpp/rcpp-and-c")
library("cRcpp")

Implementation Details

The goal of this project is to show how a hand-written C routine (one that talks to R through R’s C API and the SEXP type) can be called from within a C++ function that is, in turn, exposed to R with Rcpp. The C code never has to know that Rcpp exists; the C++ layer is responsible for translating R objects in and out of the C routine. For comparison, the package also ships a pure Rcpp implementation of the same algorithm so the two approaches can be placed side by side.

The example algorithm is a discrete convolution of two numeric vectors a and b, which produces a vector of length length(a) + length(b) - 1.

.
├── DESCRIPTION                         # Package metadata
├── NAMESPACE                           # Function and dependency registration
├── R                                   # R functions
│   ├── cRcpp-package.R                 # Package documentation
│   └── RcppExports.R                   # Autogenerated R to C++ bindings by Rcpp
├── README.md
├── man                                 # Package documentation
│   ├── convolve_cpp.Rd
│   ├── convolve_from_c.Rd
│   └── cRcpp-package.Rd
└── src                                 # Compiled code
    ├── convolve_in_c.c                 # The C routine that uses R's C API
    ├── convolve_in_c.h                 # Header exposing the C routine's prototype
    ├── convolve_from_c_to_rcpp.cpp     # Rcpp wrapper that calls the C routine
    ├── convolve_in_cpp.cpp             # Equivalent algorithm written purely in Rcpp
    └── RcppExports.cpp                 # Autogenerated R bindings

The C Routine

The C routine lives in src/convolve_in_c.c and is written directly against R’s C API. It receives its arguments as SEXP objects, coerces them to real (double) vectors, protects them from garbage collection with PROTECT(), allocates an output vector of length na + nb - 1, and fills it with the convolution. The matching UNPROTECT() releases the protected objects before the result is returned.

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

#include "convolve_in_c.h"

SEXP convolve_c(SEXP a, SEXP b) {
  int na, nb, nab;
  double *xa, *xb, *xab;
  SEXP ab;
  a = PROTECT(coerceVector(a, REALSXP));
  b = PROTECT(coerceVector(b, REALSXP));
  na = length(a); nb = length(b);
  nab = na + nb - 1;
  ab = PROTECT(allocVector(REALSXP, nab));
  xa = REAL(a); xb = REAL(b); xab = REAL(ab);
  for(int i = 0; i < nab; i++)
    xab[i] = 0.0;
  for(int i = 0; i < na; i++)
    for(int j = 0; j < nb; j++)
      xab[i + j] += xa[i] * xb[j];
  UNPROTECT(3);
  return ab;
}

The routine’s prototype is declared in a small header, src/convolve_in_c.h, so that both the C source above and the C++ wrapper below can see the same declaration.

#ifndef CONVOLVE_C_H
#define CONVOLVE_C_H

SEXP convolve_c(SEXP a, SEXP b);

#endif /* CONVOLVE_C_H */

Calling C from Rcpp

The bridge between the C routine and R is supplied in src/convolve_from_c_to_rcpp.cpp. Because the routine is compiled as C, its header must be included inside an extern "C" block so the C++ compiler uses C linkage and does not mangle the convolve_c symbol. The exported function convolve_from_c() passes its Rcpp::NumericVector arguments straight to convolve_c() (they convert implicitly to SEXP) and then wraps the returned SEXP back into an Rcpp::NumericVector. The // [[Rcpp::export]] attribute is what tells Rcpp to generate the R-level binding.

#include "Rcpp.h"

// Define the method signature

#ifdef __cplusplus
extern "C" {
#endif

#include "convolve_in_c.h"

#ifdef __cplusplus
}
#endif

// [[Rcpp::export]]
Rcpp::NumericVector convolve_from_c(const Rcpp::NumericVector& a,
                                    const Rcpp::NumericVector& b) {

  // Compute the result in _C_
  // Import it into R
  SEXP ab = convolve_c(a, b);

  // Cast as a NumericVector
  Rcpp::NumericVector result( ab );

  // Return result
  return result;
}

The exported C++ functions also carry //' roxygen comments (trimmed from the excerpts here for brevity); Rcpp::compileAttributes() copies these into R/RcppExports.R, which is where the package’s R help pages are generated from. From R the function is then called directly:

convolve_from_c(1:5, 5:1)

A Pure Rcpp Equivalent

For comparison, src/convolve_in_cpp.cpp implements the same convolution entirely with Rcpp types, without dropping down to R’s C API. This is the more idiomatic Rcpp approach and is useful for contrasting the bookkeeping required by the C routine (manual PROTECT()/UNPROTECT(), SEXP coercion) against the convenience of Rcpp::NumericVector.

#include "Rcpp.h"

// [[Rcpp::export]]
Rcpp::NumericVector convolve_cpp(const Rcpp::NumericVector& a,
                                 const Rcpp::NumericVector& b) {

  // Declare loop counters, and vector sizes
  int i, j,
  na = a.size(), nb = b.size(),
  nab = na + nb - 1;

  // Create vector filled with 0
  Rcpp::NumericVector ab(nab);

  // Crux of the algorithm
  for(i = 0; i < na; i++) {
    for(j = 0; j < nb; j++) {
      ab[i + j] += a[i] * b[j];
    }
  }

  // Return result
  return ab;
}
convolve_cpp(1:5, 5:1)

DESCRIPTION

Surfacing the C++ code with Rcpp requires Rcpp to be listed under both LinkingTo (for the headers used at compile time) and Imports (so it is available at run time).

LinkingTo:
    Rcpp
Imports:
    Rcpp (>= 1.0.12)

Source

The full package source is available on GitHub, and the rendered package site is at rd-rcpp.thecoatlessprofessor.com/rcpp-and-c.