<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-2-c-vs-fortran-vs-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/23/the-need-for-speed-part-2-c-vs-fortran-vs-c/"> 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>In my <a rel="nofollow" href="https://www.avrahamadler.com/2018/12/09/the-need-for-speed-part-1-building-an-r-package-with-fortran/" target="_blank">previous post</a>, I described the method I use for compiling Fortran (or C) into an R package using the <code>.Call</code> interface. This post will compare the speed of various implementations of the layer loss cost function.</p>
<h1>The Function</h1>
<p>Often, insurance or reinsurance is bought in stratified horizontal layers. For example, an auto policy with a $300,000 limit may have its middle third insured. This means that the primary insurer would pay up to the first $100,000 of loss, a reinsurer would pay the next $100,000, and the primary insurer would pay the last $100,000. In actuarial terms, we say that the reinsurer is covering $100K excess of $100K. The reinsured losses attach at $100K and are limited to a total of $100K paid loss. If <img src="https://i0.wp.com/www.avrahamadler.com/wp-content/ql-cache/quicklatex.com-497c874112e74efdba278710c5e1892b_l3.png?resize=18%2C14&ssl=1" alt="X" width="18" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> is the random variable representing the loss, <img src="https://i1.wp.com/www.avrahamadler.com/wp-content/ql-cache/quicklatex.com-f8016ff830b491e0b1f3122a41ccff3f_l3.png?resize=14%2C14&ssl=1" alt="L" width="14" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> is the limit of the policy, and <img src="https://i1.wp.com/www.avrahamadler.com/wp-content/ql-cache/quicklatex.com-3fae06210cb4529155c990d8d9838d68_l3.png?resize=16%2C14&ssl=1" alt="A" width="16" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="14"> the attachment, then the layer loss cost, LLC is:</p>
<p><span>   </span><span>   </span><img src="https://i0.wp.com/www.avrahamadler.com/wp-content/ql-cache/quicklatex.com-ec4441d4c42869ce26c6f3f01cef5c64_l3.png?resize=283%2C23&ssl=1" width="283" alt="[ LLC = max(0, min(X - A, L)) ]" title="Rendered by QuickLaTeX.com" data-recalc-dims="1" height="23"></p>
<p>In R, this is extremely easy to write:</p>
<pre>LLC_r &lt;- function(x, l, a) sum(pmax(0, pmin(x - a, l)))</pre>

<h1>The Question</h1>
<p>The two main languages which comprise the compiled code in R are C/C++ and Fortran. I’m not sure why, but I have a soft spot in my heart for Fortran, specifically the free-form, actually readable, versions of Fortran from F90 and on. For the packages I maintain on CRAN, I’ve found that, more often than not, Fortran will outperform C++ by a small margin—although not always. This was opportunity to directly compare the performance of C++ and Fortran.</p>
<h2>Parallelization</h2>
<h3>Basics</h3>
<p>Speed improvements often come from explicit parallelization of the compiled code. In parallel computing, <strong>vectorization</strong> means using the multiple SIMD lanes in modern chips (SSE/AVX/AVX2). The loop is only performed in one thread but it can be chunked inside the CPU. Consider a one-lane highway approaching a toll plaza; the speed limit remains the same but just at the tollbooth four cars can be processed at the same time. <strong>Parallelization</strong> means taking advantage of multiple cores. This is akin to splitting the one lane into several lanes, each one with its own speed limit and approach to its own toll-plaza. Sometimes single-thread vectorization is called “fine-grain parallelism” and multi-thread parallelism is called “coarse-grain parallelism”. Both can be combined on today’s modern SIMD CPUs with multiple cores.</p>
<p>In theory, either of these should be faster than single-thread serial calculation. In practice, however, there is a certain amount of overhead involved in setting up either kind of parallelism. Therefore, for small loads, it is not uncommon to find the serial versions to be the fastest. Even for larger, or huge, loads, it is not always the case that the most parallelism is worth the most overhead. In each language tested, the code will test coarse-grained, fine-grained, and combined parallelism.</p>
<h3>R-specific hurdle</h3>
<p>The standard way to parallelize C/C++/Fortran is by using <a rel="nofollow" href="https://www.openmp.org/" target="_blank">OpenMP</a>. However, R’s mechanics officially don’t allow OpenMP to be used by multiple languages in the same package. While it may be forced, the forcing violates R norms and relies on some default linker flags being the same across languages. Even now, December 2018, the Fortran mechanics for OpenMP are changing. So to test the languages, I wrote two small packages—one exclusively C++ using Rcpp and one using the C/Fortran interface. In the end, I did violate those norms to show some pure C calls which were put into the Fortran package. While this will work for now and for this post, it is not good practice and should be avoided.</p>
<h1>Approaches</h1>
<h2>R</h2>
<p>In base R, the pure translation would be the nested use of <code>pmin</code> and <code>pmax</code>. There is one variant; the documentation for the two functions also describes a <code>.int</code> version for both which are “faster internal versions only used when all arguments are atomic vectors and there are no classes.” So for R, there are two options to be checked:</p>
<pre>LLC_r1 &lt;- function(x, l, a) sum(pmax(0, pmin(x - a, l)))
LLC_r2 &lt;- function(x, l, a) sum(pmax.int(0, pmin.int(x - a, l)))</pre>

<h2>C++/Rcpp</h2>
<h3>Serial Loop</h3>
<p>Unlike R, the easiest way to calculate this function in C++ is to use a for loop. There may be more elegant ways using STL, but those are more difficult to implement (and the LLC function isn’t unary anyway). In R, loops are deprecated in favor of vectorized functions and the <code>apply</code> family. The reason is that the looping is handled by the underlying C/C++/Fortran and not at the R level. Here we are living in the C++ basement, so well-formed loops are both good and fast. The archetypal loop implementation for the layer loss cost function in Rcpp would be:</p>
<pre>// [[Rcpp::export]]
double LLC_cpp(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>Not only is this implementation fast, a sufficiently good and aggressive compiler can also attempt autovectorization for even more speed (<code>-O3</code> will usually attempt this). However, explicit vectorization and parallelization for C/C++ and Fortran will be applied using OpenMP. There are other ways to impose parallelism, like the Intel TBB framework behind RcppParallel. In that lies a possible addendum or two to this post.</p>
<h3>Parallel Loops</h3>
<p>Three additional flavors of the C++ code were tested:</p>
<ul>
<li>Only using SIMD vectorization: <code>LLC_cpp<strong>s</strong></code></li>
<li>Only using explicit parallization: <code>LLC_cpp<strong>l</strong></code></li>
<li>Using both: <code>LLC_cpp<strong>ls</strong></code></li>
</ul>
<p>The entire <code>.cpp</code> file is posted below to show the includes and other Rcpp necessities which weren’t covered in the first post.</p>
<pre>#include &lt;Rcpp.h>
#include </pre>