The stochastic logistic population model described in this blog post has become my default exercise when I'm exploring a new programming language (Racket, F#, Rust). I ended that Rust post by noting that I was interested in Rust as a modern alternative for scientific computing and thought it would be a good learning exercise to re-write small legacy Fortran programs in Rust. In the process of looking for Fortran programs to translate to Rust, though, I found myself becoming more interested in the idea of learning Fortran than Rust, particularly after learning about the efforts to improve the tooling around Fortran (e.g., here and here). So, here we are...exploring Fortran via the stochastic population model exercise.
Standalone Fortran Program
First, I wrote a standalone Fortran program of the stochastic population model. The main thing that tripped me up was reading arguments from the command line. I initialized a 1D array for holding the values of the arguments.
character(len=20), dimension(5) :: args
I initially missed the
len argument for
character and was confused by the results I was getting until I realized that the program was only reading the first character provided for each argument.
We can straightforwardly fill the
args array with a loop.
do i=1,command_argument_count() call get_command_argument(i, args(i)) end do
But the syntax for converting the arguments from strings to the correct types is unusual.
read(args(1), *) t
t is the first value in
args (indexed with
* indicates the default format of the value, which the compiler infers is an
integer because that is the type of
I also created a file that only contains subroutines for calling from R. Two of the subroutines (
rnormal) are mostly the same as from the main program, but the names are changed because the R help for the
.Fortran function recommended not using underscores in the subroutine names and the type of several parameters is
double precision, not
real, for compatibility with R. One of the subroutines,
rnorm, is only included as a test of
subroutine rnorm(reps, mu, sigma, arr) implicit none integer :: i integer, intent(in) :: reps double precision, intent(in) :: mu, sigma double precision, intent(inout) :: arr(reps) do i=1,reps call rnormal(mu, sigma, arr(i)) end do end subroutine rnorm
arr parameter is declared as
intent(inout) because on the R side we will pass in a vector that is filled on the Fortran side and then returned again to R. There are performance pitfalls with this approach if you are passing large vectors back-and-forth between Fortran and R, but those can be partly mitigated with the use of the dotCall64 package.
We compile this file into a shared object for use in R with
R CMD SHLIB stochastic-logistic-subroutines.f90
Calling Fortran from R
In the R code, we load the shared object with
To test that the
rnormal subroutine is giving reasonable results, we call the
rnorm subroutine as the first argument to
.Fortran. In this example, we draw 1 million random variates from a normal distribution with a mean of 7.8 and standard deviation of 8.3.
> n = 1000000L > rnorm_result <- .Fortran("rnorm", n, 7.8, 8.3, vector("numeric", n)) # .Fortran returns a list; 4th element of list in this case contains the vector > mean(rnorm_result[])  7.798748 > sd(rnorm_result[])  8.302705
Following the original blog post, we benchmark the R, Rcpp, and Fortran versions.
> library(microbenchmark) > microbenchmark( Rcpp = logmodc(t, yinit, r, k, thetasd), Fortran = .Fortran("logmodf", t, yinit, r, k, thetasd, vector("numeric", t)), R = logmodr(t, yinit, r, k, thetasd), times = 500L ) Unit: microseconds expr min lq mean median uq max neval Rcpp 13.235 14.6260 17.11903 15.5280 16.8100 86.612 500 Fortran 21.950 22.8965 26.76053 23.8760 25.3160 72.528 500 R 56.994 58.8775 67.25860 59.9085 62.3655 243.965 500
The original blog post (from 2017) observed a 19x speed up of the Rcpp version over the R version. Here we only get a 4x improvement presumably related to performance improvements in R. Our Fortran version is 2.5x faster than the R version.
Again, following the original blog post (with a slight modification), we benchmark running multiple simulations.
> rseq <- seq(1.1, 2.2, length.out = 10) > microbenchmark( Rcpp = lapply(rseq, function(x) logmodc(t, yinit, x, k, thetasd)), Fortran = lapply(rseq, function(x) .Fortran("logmodf", t, yinit, x, k, thetasd, vector("numeric", t))), R = lapply(rseq, function(x) logmodr(t, yinit, x, k, thetasd)), times = 500L ) Unit: microseconds expr min lq mean median uq max neval Rcpp 159.578 166.4295 214.4782 198.2240 202.9475 8265.555 500 Fortran 233.397 239.9620 264.1335 256.8515 262.1665 2994.857 500 R 621.029 635.1115 676.0114 665.7525 674.2735 2952.635 500
The relative performance gains for multiple simulations are similar to the single simulation.
This post doesn't make a resounding case for using Fortran over C++ for speeding up simulations in R. The C++ code is shorter and runs faster. My C++ avoidance has centered on the complexity of the language. Perhaps there is a simpler subset of C++ that would go a long way towards speeding up my R code. However, this post makes the case that Fortran is still easier to learn than C++ for physics simulations.