<p><div><div><div><div data-show-faces="false" data-height="21" data-width="100" data-layout="button_count" data-send="true" data-href="https://www.r-bloggers.com/the-need-for-speed-part-1-building-an-r-package-with-fortran-or-c/"></div></div><div></div></div><div></div></div>
<div>
(This article was first published on <strong><a href="https://www.avrahamadler.com/2018/12/09/the-need-for-speed-part-1-building-an-r-package-with-fortran/"> R – Strange Attractors</a></strong>, and kindly contributed to <a rel="nofollow" href="https://www.r-bloggers.com/">R-bloggers)</a> </div>
<p><a rel="nofollow" href="https://commons.wikimedia.org/wiki/File:Algol%26Fortran_family-by-Borkowski.svg" target="_blank"><img src="https://i2.wp.com/upload.wikimedia.org/wikipedia/commons/thumb/8/84/Algol%26Fortran_family-by-Borkowski.svg/288px-Algol%26Fortran_family-by-Borkowski.svg.png?resize=288%2C240&ssl=1" alt="Genealogy tree of programming languages - Algol &amp; Fortran family " width="288" data-recalc-dims="1" height="240"></a></p>
<h1>Searching for Speed</h1>
<p>Everyone who has ever used R has, at one time or another, wished for an increase in R’s speed. If you haven’t, you’re not using R hard enough! Recently, as part of some research on credibility, I was calculating layer loss costs for millions of simulated loss observations. As I progressed, the R markdown document I was was creating took longer and longer to render. Profiling the code showed one of the bottlenecks was the layer loss cost calculation. Being actively interested in compiled code for R for the better part of a decade, the simplicity of the function was a great opportunity for me to compare different implementations of compiled code.</p>
<p>I started writing one long post, but realized it would be better split into two parts. This post will describe the important steps in building an R package using Fortran. The next post will describe and compare the form and timings of the base R, C++, and Fortran implementations of the layer loss cost function. I hope that compiling (bad pun most definitely intended) this information and guidance in one place will make it easier for people to develop their own enhancements to R. Remember, a package doesn’t have to mean CRAN-ready code—a package is a collection of code that is easy to load and use in R. Don’t be afraid to experiment and push your own R boundaries!</p>
<h1>Why Fortran?</h1>
<p>In my opinion, Fortran gets a bad rap. It remains the oldest compiled language in active use. Yes, you wouldn’t want to write a web browser or a CRPG in it, but for mathematical calculation it remains difficult to beat <strong>For</strong>mula <strong>tran</strong>slator as a language. The majority of high-performance computing today remains in Fortran—it has been proven to work for the past 60+ years, it works well, and free-form Fortran is relatively easy to read. Searching for “why is Fortran still relevant” will <a rel="nofollow" href="http://moreisdifferent.com/2015/07/16/why-physicsts-still-use-fortran/" target="_blank">return</a> <a rel="nofollow" href="https://www.informationweek.com/software/enterprise-applications/fortran-7-reasons-why-its-not-dead/d/d-id/1321174" target="_blank">many</a> <a rel="nofollow" href="https://hackernoon.com/fortran-is-still-a-thing-fc84df4cf638" target="_blank">answers</a>, even in today’s world of Rust, Python, and Julia.</p>
<h1>Fortran in R</h1>
<p>R itself uses Fortran for much of its hard number crunching. There is a built-in base-R interface named <code>.Fortran</code>. However, according to <a rel="nofollow" href="https://github.com/wrathematics" target="_blank">Drew Schmidt</a>, an <a rel="nofollow" href="http://librestats.com/hpcdevcon2016/files/tutorial/hpcR.pdf" target="_blank">expert</a> in <a rel="nofollow" href="https://www.nics.tennessee.edu/r-tutorial-worldwide-2015" target="_blank">high-performance R</a> computation, this interface is outdated, comes with significant overhead, and should not be used if at all possible. His comments on his <a rel="nofollow" href="https://github.com/wrathematics/Romp/blob/master/README.md" target="_blank">Romp package on github</a>, and the associated files, are very valuable in learning how to navigate the C/Fortran API for R. My <a rel="nofollow" href="https://cran.r-project.org/package=Delaporte" target="_blank">Delaporte</a> package’s conversion from C++ to C/Fortran, and the resulting speedup, was heavily influenced by his writing.</p>
<p>There are <a rel="nofollow" href="http://users.stat.umn.edu/~geyer/rc/" target="_blank">other</a> <a rel="nofollow" href="https://rollingyours.wordpress.com/2014/04/11/fortran-and-r-speed-things-up/" target="_blank">pages</a> discussing the <code>.Fortran</code> interface. Therefore, this post will focus on how to use use Fortran by building it into a package and using <code>.Call</code>. Readers will benefit more if they have basic knowledge of R package generation (NAMESPACE, DESCRIPTION, /man, etc.). For those who do not, reading the links below should prove valuable for basic package generation as well.</p>
<h1>Before writing any code</h1>
<h2>Necessary pre-reading</h2>
<p>At first glance, writing compiled code in R packages using the <code>.Call</code> interface is complicated. To say that the R API is obtuse and not-well-described would be charitable. It can be downright byzantine. This is why <a rel="nofollow" href="http://dirk.eddelbuettel.com/" target="_blank">Dirk Eddelbuettel</a>’s Rcpp package changed the world. It allowed authors to focus on the calculation code and leave the interfacing to Rcpp. As there is no Rcpp-like package for pure C or Fortran, some basic knowledge of the C API is necessary. The following links are not just highly recommended, they should be considered required reading:</p>
<ul>
<li><a rel="nofollow" href="https://cran.r-project.org/doc/manuals/r-release/R-exts.html" target="_blank">Writing R Extensions</a>: Read this a minimum of six times from start to finish. Re-read it before every time you plan to ask a question of <a rel="nofollow" href="https://stat.ethz.ch/mailman/listinfo/r-package-devel" target="_blank">R-package-devel</a>.</li>
<li><a rel="nofollow" href="http://www.stat.purdue.edu/~liu105/STAT598G_lab/Rcourse94up.pdf" target="_blank">Dr. Ripley’s slides</a> on calling other languages from R: Over a decade old now, but the mechanics are the same.</li>
<li><a rel="nofollow" href="http://adv-r.had.co.nz/C-interface.html" target="_blank">Hadley Wickham on the C API</a>: Hadley has forgotten more about R than most of us will ever know.</li>
<li><a rel="nofollow" href="https://tolstoy.newcastle.edu.au/R/e17/devel/att-0724/R_API_cheat_sheet.pdf" target="_blank">Simon Urbanek’s cheat sheet</a>: Dense, but chock full of information. Once you have some idea, this is a good reminder/reference.</li>
<li><a rel="nofollow" href="http://mazamascience.com/WorkingWithData/?p=1099" target="_blank">Jonathan Callahan</a> on <code>.C</code> vs. <code>.Call</code>: Applies to <code>.Fortran</code> as well. See Simon Urbanek’s response to Hervé Pagès at the beginning.</li>
<li><a rel="nofollow" href="https://raw-r.org/R_API.php" target="_blank">Ben Bass on the C API</a>: This wasn’t around when I started learning about the API; I wish it had been.</li>
</ul>
<h2>Critical points to remember</h2>
<ul>
<li>Use the <code>.Call</code> interface. It’s modern, has lower overhead, and is safer than the <code>.C</code>, <code>.Fortran</code>, or <code>.External</code> interfaces.</li>
<li>Everything that moves between R and C is a SEXP. <em><strong>Everything</strong></em>. If it stays in C, C++, or Fortran, or moves between those, that’s OK, but R only understands SEXPs. It’s fine to not understand the deep mechanics of an <strong>S EXP</strong>pression pointer. Just know that data moving between R and C via <code>.Call</code> must be defined as SEXPs in C. Inside the function is where the actual type will be specified.</li>
<li>After GCC 4.5, it’s easy to call Fortran from C but only as a Fortran <strong>subroutine</strong> and not a function. Inside the Fortran code functions are fine, but any Fortran code that is going to be called by C and then passed to R should be a subroutine.</li>
<li>Read “Writing R Extensions” another five or ten times.</li>
</ul>
<h1>How to build a package using compiled code (Fortran)</h1>
<h2>Brief Description of the Function</h2>
<p>The function precipitating these posts is the insured layer loss cost. This is a simple function whose implementation in base R is one line:</p>
<pre>LLC_r &lt;- function(x, l, a) sum(pmax(0, pmin(x - a, l)))</pre>

<p>In C++ (Rcpp-style), the simplest implementation would be a loop:</p>
<pre>// [[Rcpp::export]]
double LLC_c(NumericVector X, double L, double A) {
  int n = X.size();
  double LLC = 0.0;
  for (int i = 0; i &lt; n; ++i) {
    LLC += fmax(0.0, fmin(X[i] - A, L));
  }
  return(LLC);
}</pre>

<p>You will have to include <code>cmath</code> to get <code>fmax</code> and <code>fmin</code>. We’ll see this again in the next post.</p>
<h2>First step: Write the Fortran</h2>
<p>I’ve found it easiest to create Fortran files as modules. This allows simpler calling of functions and subroutines from different files, although it does require some care in writing the Makevars if there will be dependencies between separate Fortran files. For the actual calculation, the direct analogue of the loop in C++ would be a loop in Fortran. For this example, the complete Fortran file with some added comments looks like this:</p>
<pre>module fortloop
    use, intrinsic :: iso_c_binding

    implicit none
    private
    public :: llc_f

contains

    subroutine llc_f(x, n, l, a, llc) bind(C, name = &quot;llc_f_&quot;)

        real(kind = c_double), intent(in)               :: l, a !limit &amp; attach
        integer(kind = c_int), intent(in), value        :: n    !Length of x
        real(kind = c_double), intent(in), dimension(n) :: x    !Vector of loss
        real(kind = c_double), intent(out)              :: llc  !Output variable
        integer                                         :: i    !Internal count

        llc = 0.0_c_double
        do i = 1, n
            llc = llc + max(0.0_c_double, min(x(i) - a, l))
        end do

    end subroutine llc_f

end module fortloop</pre>

<p>The subroutine is passed a vector <code>x</code> of length <code>n</code>, the limit <code>l</code>, the attachment <code>a</code>, and the return variable <code>llc</code>. Using <code>iso_c_binding</code> allows the coder to explicitly state the name of the calling C program in the subroutine via <code>bind</code>. The <code>F77_</code> prefixes in C add a trailing underscore, so remember to add it to the bound name. <code>iso_c_binding</code> also has built-in variable-type definitions; use them for easier interoperability. Using <code>implicit none</code>—thus needing to define the variable types—is helpful in ensuring well-formed code. Explicitly stating which routines are public and private is a matter of preference.</p>
<p>Save this file into the package’s <code>src</code> directory as <code>./src/fortloopf.f95</code>. Why <code>.f95</code>? If you cannot answer that question, apparently you haven’t read Writing R Extensions enough! The next step is to write the C function which will call Fortran and the seperate C function which will communicate with R.</p>
<h2>Second step: Write the C</h2>
<h3>C-to-Fortran</h3>
<p>The first function is the C code which will call Fortran, and it’s rather simple:</p>
<pre>void F77_NAME(llc_f)(double *x, int n, double *l, double *a, double *ret);</pre>

<p>This function passes the variables to Fortran in the order it expects. This was the function “bound” in Fortran with a trailing underscore. Almost always, variables are passed to Fortran by reference, and thus need to be pointers in C. Fortran 2003 introduced the <code>value</code> keyword, which allows passing by value. In the Fortran subroutine written above, <code>n</code> was declared with the <code>value</code> keyword, so in C it does not need to be a pointer. The <code>F77_NAME</code> prefix is explained in <a rel="nofollow" href="https://cran.r-project.org/doc/manuals/r-release/R-exts.html#Calling-C-from-FORTRAN-and-vice-versa" target="_blank">Writing R extensions</a>, which you should have read near a dozen times by now!</p>
<h3>C-to-R</h3>
<p>The second function is the one which will be called from R using the <code>.Call</code> interface:</p>
<pre>extern SEXP c_llc_f(SEXP x, SEXP l, SEXP a){
  const int n = LENGTH(x);
  SEXP ret;
  PROTECT(ret = allocVector(REALSXP, 1));
  F77_CALL(llc_f)(REAL(x), n, REAL(l), REAL(a), REAL(ret));
  UNPROTECT(1);
  return(ret);
}</pre>

<p>This function communicates between R and C. Therefore, the vector of losses <code>x</code>, the limit <code>l</code>, and the attachment <code>a</code> need to be SEXPs. The calculation of the length of <code>x</code> is done in C via <code>LENGTH</code> and assigned to <code>n</code>. Since it isn’t coming from R, it doesn’t need to be a SEXP. The return value, <code>ret</code>, <strong>is</strong> going back to R, so it needs to be defined as a SEXP. Since it is will be a double precision real variable, its flavor of SEXP is REALSXP. The <code>PROTECT</code> and <code>UNPROTECT</code> calls are needed to ensure that R’s garbage collection doesn’t free the memory until we are done with it. Since the Fortran loop returns a single value, <code>ret</code>(urn) is declared as SEXP and given protected allocation of a length 1 real vector. Remember, there are no scalars in R.</p>
<p>In the actual call (<code>F77_CALL</code> calling the <code>F77_NAME</code>), the SEXPs have to be appropriately typed. Since they are all doubles, they get cast via the <code>REAL()</code> call. The length variable, <code>n</code>, is not being passed in from R. Since it was defined inside the function as a <code>const int</code>, it doesn’t need casting in the call to the C-to-Fortran function. If it was an SEXP that came from R, it <em>would</em> need casting as <code>INTEGER(n)</code>. In this case it would also need to be a pointer to that value as well. There are other types of SEXP, of which I’m sure you are now aware, having read all the links above.</p>
<p>Save these two functions into the package’s <code>src</code> directory as <code>./src/fortloopc.c</code>. This file is incomplete, but the final C coding will be easier to understand once the R code is written.</p>
<h3>Hard-learned lessons</h3>
<h4>Booleans</h4>
<p>One valuable point, gleaned from <a rel="nofollow" href="https://stat.ethz.ch/pipermail/r-package-devel/2017q1/001383.html" target="_blank">much frustration</a>, relates to passing logicals. <code>iso_c_binding</code> includes boolean bindings. For almost all Fortran compilers, this would allow easy passing of a logical out of R into a boolean in C, followed by the passing of that boolean into a logical Fortran variable. Sadly, that does not work with some Fortran implementations—specifically Solaris SPARC which has opposite endianness from most other CPUs. The foolproof way to pass booleans from R to Fortran via C is to cast them as integers—0 for FALSE and 1 for TRUE is the usual convention—and then test for equality to 0 or 1 in Fortran. Maybe one day all major implementations will work seamlessly. Until then, while it is less elegant, passing integers can save a whole host of <a rel="nofollow" href="https://www.avrahamadler.com/2017/03/31/delaporte-package-the-sparcmonster-is-sated/" target="_blank">headaches</a>.</p>
<h4>Portability</h4>
<p>If your intent is to eventually create a CRAN-worthy package, pay special attention to portability. Make sure you know which elements of Fortran 90, 95, 2003, or 2008 are implemented in the regular compilers. For example, as of December 2018, the oldest CRAN compiler is <a rel="nofollow" target="_blank">GCC 4.6.3 for r-oldrel-windows</a>. So to be compatible with <code>r-oldrel-windows-ix86+x86_64</code> Fortran cannot use the intrinsic function<code>log_gamma</code>, which was only added in Fortran 2008. The user must write their own log gamma code or use an existing module.</p>
<p>Another example is that the current release of R on Windows can use Fortran with <code>iso_c_binding</code> for variable types and subroutines. That was <a rel="nofollow" href="https://gcc.gnu.org/wiki/GFortran/News#GCC4.3" target="_blank">implemented in GCC 4.3</a> and is in Oracle Developer Studio 12.5 for Solaris. However, the current release of R on Windows cannot use <code>ieee_arithmetic</code>, as that was only implemented in <a rel="nofollow" href="https://gcc.gnu.org/wiki/GFortran/News#gfortran_5" target="_blank">GFortran 5</a>. The current Windows toolchain is based on version 4.9.3, although there is work on updating the Windows toolchain to GCC 8. </p>
<h2>Third step: Write the R</h2>
<p>In R, the function to call the C code is:</p>
<pre>LLC_f2 &lt;- function(x, l, a) {
  if (!is.double(x)) {storage.mode(x) &lt;- 'double'}
  if (!is.double(l)) {storage.mode(l) &lt;- 'double'}
  if (!is.double(a)) {storage.mode(a) &lt;- 'double'}
  .Call(c_llc_f, x, l, a)
}</pre>

<p>The user calls this function the same way they would the base R version. The mechanics are to first ensure that the variables are of the right type (passing the wrong type from R to C has been known to release a horde of demons somewhere within a ten mile radius of the infraction) and then <code>.Call</code>ing the C. Save this under the package’s <code>R</code> directory as <code>./R/LLC.R</code>.</p>
<h2>Fourth step: Add the needed package code to C</h2>
<p>Recently, R-core has become stricter with registration of native routines for compiled code. But you knew that because you follow <a rel="nofollow" href="https://stat.ethz.ch/mailman/listinfo/r-devel" target="_blank">R-devel</a> and have all-but memorized Writing R Extensions at this point, right? So you already know that there are two steps: the first is to declare the interface methods and the second is to register the routines. In this case, there is only use of the <code>.Call</code>, so the registration will only require one type. Each function needs to be defined with the number of variables being passed between R and C.</p>
<pre>static const R_CallMethodDef CallEntries[] = {
  {&quot;c_llc_f&quot;,   (DL_FUNC) &amp;c_llc_f,   3},
  {NULL,         NULL,                0}
};

void R_init_SimpFort (DllInfo *dll) {
  R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
  R_useDynamicSymbols(dll, FALSE);
}</pre>

<p>The first function creates the definitions for all functions using the <code>.Call</code> interface—the only interface you should ever use—and the second registers the routines and passes FALSE to Dynamic Symbols, which is highly preferred by R. The entire C code file, <code>fortloop.c</code> now looks like this:</p>
<pre>#include &lt;R.h>
#include &lt;Rinternals.h>
#include &lt;stdlib.h> // for NULL
#include &lt;Rmath.h>
#include </pre>