Extending R with C++ and Fortran
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
A recent social-media
question by James
Curran inquired about the best, or recommended ways, to
extend R with Fortran code. Part of the question was whether the .Fortran()
interface was still
recommended or as there is ‘conflicting advice’ out there. Dirk
then followed up and pointed to the
(stunning!) performance gains reported by
glmnet
which switched from
.Fortran()
to a C++ interface using Rcpp and the (now much preferred) .Call()
interface. One
key reason behind the performance gains is that .Fortran()
requires copies of all arguments, just
like the (also effectively deprecated) .C()
interface. Whereas .Call()
works with SEXP
object
which are pointers: this can be dramatically faster and more efficient as object sizes increase.
A few years earlier, and for a related question, JBrandon Duck-Mayr had written a very comprehensive answer on StackOverflow. It is backed by an example package mixedlang which implements the recommendation.
It starts from a Fortran90 function multiplying two ‘real’ aka double
valued inputs:
REAL*8 FUNCTION MULTIPLY (X, Y) REAL*8 X, Y MULTIPLY = X * Y RETURN END
This can be connected quite easily to C++ code using the common extern "C"
declaration (specifying
that a C calling convention is used from the C++ code). It still shows the Rcpp::depends()
used
when sourceCpp()
-ing a function, it is not needed in a package like mixedlang
.
#include "RcppArmadillo.h" // [[Rcpp::depends(RcppArmadillo)]] // First we'll declare the MULTIPLY Fortran function // as multiply_ in an extern "C" linkage specification // making sure to have the arguments passed as pointers. extern "C" { double multiply_(double *x, double *y); } // Now our C++ function // [[Rcpp::export]] Rcpp::NumericVector test_function(Rcpp::NumericVector x) { // Get the size of the vector int n = x.size(); // Create a new vector for our result Rcpp::NumericVector result(n); for ( int i = 0; i < n; ++i ) { // And for each element of the vector, // store as doubles the element and the index double starting_value = x[i], multiplier = (double)i; // Now we can call the Fortran function, // being sure to pass the address of the variables result[i] = multiply_(&starting_value, &multiplier); } return result; }
Once both functions are compiled and loaded (as e.g. in package mixedlang
) the wrapper function
can be called from R as usual:
mixedlang::test_function(0:9) # [1] 0 1 4 9 16 25 36 49 64 81
We hope the (recently updated) package at GitHub serves as starting point for other wanting to combine R and Fortran via Rcpp.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.