|
| 1 | +--- |
| 2 | +title: "Extending R with C++ and Fortran" |
| 3 | +author: "Dirk Eddelbuettel and JBrandon Duck-Mayr" |
| 4 | +license: GPL (>= 2) |
| 5 | +tags: basics |
| 6 | +summary: Rcpp can also be helpful in providing glue code to connect to Fortran. |
| 7 | +layout: post |
| 8 | +src: 2023-02-05-Combining-R-with-Cpp-and-Fortran.Rmd |
| 9 | +--- |
| 10 | + |
| 11 | +A recent social-media |
| 12 | +[question ](https://mastodon.social/@[email protected]/109770574777963511) by [James |
| 13 | +Curran](https://profiles.auckland.ac.nz/j-curran) inquired about the best, or recommended ways, to |
| 14 | +extend R with Fortran code. Part of the question was whether the `.Fortran()` interface was still |
| 15 | +recommended or as there is 'conflicting advice' out there. [Dirk](https://dirk.eddelbuettel.com) |
| 16 | +then [followed up](https://mastodon.social/@eddelbuettel/109772801398410192) and pointed to the |
| 17 | +(stunning!) performance gains [reported by |
| 18 | +`glmnet`](https://cran.r-project.org/web/packages/glmnet/news/news.html) which switched from |
| 19 | +`.Fortran()` to a C++ interface using Rcpp and the (now much preferred) `.Call()` interface. One |
| 20 | +key reason behind the performance gains is that `.Fortran()` requires copies of all arguments, just |
| 21 | +like the (also effectively deprecated) `.C()` interface. Whereas `.Call()` works with `SEXP` object |
| 22 | +which are _pointers_: this can be dramatically faster and more efficient as object sizes increase. |
| 23 | + |
| 24 | +A few years earlier, and for a related question, [JBrandon Duck-Mayr](https://jbduckmayr.com/) had written a _very comprehensive_ |
| 25 | +[answer on StackOverflow](https://stackoverflow.com/questions/31396802/integrate-fortran-c-with-r). |
| 26 | +It is backed by an example package [mixedlang](https://github.com/duckmayr/mixedlang) which |
| 27 | +implements the recommendation. |
| 28 | + |
| 29 | +It starts from a Fortran90 function multiplying two 'real' aka `double` valued inputs: |
| 30 | + |
| 31 | +```f |
| 32 | +REAL*8 FUNCTION MULTIPLY (X, Y) |
| 33 | +REAL*8 X, Y |
| 34 | +MULTIPLY = X * Y |
| 35 | +RETURN |
| 36 | +END |
| 37 | +``` |
| 38 | + |
| 39 | +This can be connected quite easily to C++ code using the common `extern "C"`declaration (specifying |
| 40 | +that a C calling convention is used from the C++ code). It still shows the `Rcpp::depends()` used |
| 41 | +when `sourceCpp()`-ing a function, it is not needed in a package like `mixedlang`. |
| 42 | + |
| 43 | +```cpp |
| 44 | +#include "RcppArmadillo.h" |
| 45 | + |
| 46 | +// [[Rcpp::depends(RcppArmadillo)]] |
| 47 | + |
| 48 | +// First we'll declare the MULTIPLY Fortran function |
| 49 | +// as multiply_ in an extern "C" linkage specification |
| 50 | +// making sure to have the arguments passed as pointers. |
| 51 | +extern "C" { |
| 52 | + double multiply_(double *x, double *y); |
| 53 | +} |
| 54 | + |
| 55 | +// Now our C++ function |
| 56 | +// [[Rcpp::export]] |
| 57 | +Rcpp::NumericVector test_function(Rcpp::NumericVector x) { |
| 58 | + // Get the size of the vector |
| 59 | + int n = x.size(); |
| 60 | + // Create a new vector for our result |
| 61 | + Rcpp::NumericVector result(n); |
| 62 | + for ( int i = 0; i < n; ++i ) { |
| 63 | + // And for each element of the vector, |
| 64 | + // store as doubles the element and the index |
| 65 | + double starting_value = x[i], multiplier = (double)i; |
| 66 | + // Now we can call the Fortran function, |
| 67 | + // being sure to pass the address of the variables |
| 68 | + result[i] = multiply_(&starting_value, &multiplier); |
| 69 | + } |
| 70 | + return result; |
| 71 | +} |
| 72 | +``` |
| 73 | +
|
| 74 | +Once both functions are compiled and loaded (as _e.g._ in package `mixedlang`) the wrapper function |
| 75 | +can be called from R as usual: |
| 76 | +
|
| 77 | +```r |
| 78 | +mixedlang::test_function(0:9) |
| 79 | +# [1] 0 1 4 9 16 25 36 49 64 81 |
| 80 | +``` |
| 81 | + |
| 82 | +We hope the (recently updated) package at GitHub serves as starting point for other wanting to |
| 83 | +combine R and Fortran via Rcpp. |
| 84 | + |
0 commit comments