Yutaka Masuda has a great tutorial on extending R with Fortran using R’s Foreign Function Interface (FFI).
https://masuday.github.io/fortran_tutorial/r.html
The FFI approach has two major limitations. First, it copies R objects, which hurts performance when passing a large R array to Fortran. Second, R does not natively support 64-bit integers. R will convert an integer larger than 2^31 - 1 to a double-precision real number. This feature makes it difficult to pass a large R array to Fortran because its size/shape needs to be an integer. Fortunately, we can use the R package dotCall64 to tackle these limitations.
Below is a Fortran file named test.f90 containing several functions to be exported to R. These functions are wrapped in a module and exported to R via subroutines. That is because R understands Fortran subroutines but not functions. In addition, these subroutines need to use explicit-shape arrays as their arguments because R only understands explicit-shape Fortran arrays. However, the functions in the module can use assumed-shaped arrays as their arguments. These functions are used to illustrate how to pass R arrays of various shapes to Fortran. 64-bit integer is used since you are more likely to need Fortran when working with large R arrays. OpenMP is used to showcase Fortran’s parallel computing capability.
!import function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
module test
use ieee_arithmetic, only: ieee_value, ieee_quiet_nan
private
public :: sum1d, sum2d, sum3d, lm, nan_fill
contains
function sum1d(x) result(res)
integer(8) :: n, i
real(8) :: x(:), res
n = size(x, kind = 8)
res = 0.0_8
!$omp parallel do reduction(+:res)
do i = 1, n
res = res + x(i)
end do
end function
function sum2d(x) result(res)
integer(8) :: n1, n2, i1, i2
real(8) :: x(:, :), res
n1 = size(x, dim = 1, kind = 8)
n2 = size(x, dim = 2, kind = 8)
res = 0.0_8
!$omp parallel do reduction(+:res) collapse(2)
do i2 = 1, n2
do i1 = 1, n1
res = res + x(i1, i2)
end do
end do
end function
function sum3d(x) result(res)
integer(8) :: n1, n2, n3, i1, i2, i3
real(8) :: x(:, :, :), res
n1 = size(x, dim = 1, kind = 8)
n2 = size(x, dim = 2, kind = 8)
n3 = size(x, dim = 3, kind = 8)
res = 0.0_8
!$omp parallel do reduction(+:res) collapse(3)
do i3 = 1, n3
do i2 = 1, n2
do i1 = 1, n1
res = res + x(i1, i2, i3)
end do
end do
end do
end function
function lm(x, y) result(res)
integer(8) :: n, i
real(8) :: x(:), y(:), res, sum_x, sum_y, sum_xx, sum_xy
n = min(size(x, kind = 8), size(y, kind = 8))
sum_x = 0.0_8
sum_y = 0.0_8
sum_xx = 0.0_8
sum_xy = 0.0_8
!$omp parallel do reduction(+:sum_x, sum_y, sum_xx, sum_xy)
do i = 1, n
sum_x = sum_x + x(i)
sum_y = sum_y + y(i)
sum_xx = sum_xx + x(i)*x(i)
sum_xy = sum_xy + x(i)*y(i)
end do
res = (n*sum_xy - sum_x*sum_y)/(n*sum_xx - sum_x*sum_x)
end function
function nan_fill(x, k) result(res)
integer(8) :: n, k, i
real(8) :: x(:), res(size(x, kind = 8))
n = size(x, kind = 8)
!$omp parallel do
do i = 1, k
res(i) = ieee_value(0.0_8, ieee_quiet_nan)
end do
!$omp parallel do
do i = k + 1, n
res(i) = x(i)
end do
end function
end module
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!export function
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fortran_sum1d(n, x, res)
use test, only: sum1d
integer(8) :: n
real(8) :: x(n), res
res = sum1d(x)
end subroutine
subroutine fortran_sum2d(n1, n2, x, res)
use test, only: sum2d
integer(8) :: n1, n2
real(8) :: x(n1, n2), res
res = sum2d(x)
end subroutine
subroutine fortran_sum3d(n1, n2, n3, x, res)
use test, only: sum3d
integer(8) :: n1, n2, n3
real(8) :: x(n1, n2, n3), res
res = sum3d(x)
end subroutine
subroutine fortran_lm(n, x, y, res)
use test, only: lm
integer(8) :: n
real(8) :: x(n), y(n), res
res = lm(x, y)
end subroutine
subroutine fortran_nan_fill(n, x, k, res)
use test, only: nan_fill
integer(8) :: n, k
real(8) :: x(n), res(n)
res = nan_fill(x, k)
end subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Below is an R file importing the above functions and validating the results with corresponding R codes.
#header
####################################################################################################
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
require(dotCall64)
fortran_source <- function(x)
{
#x: f90 file (character)
x_mod <- paste0(substr(x, 1, nchar(x) - 4), ".mod")
x_dynlib <- paste0(substr(x, 1, nchar(x) - 4), .Platform$dynlib.ext)
cmd <- paste("gfortran",
"-fimplicit-none",
"-O3",
"-fopenmp",
"-fpic",
"-shared",
x,
"-o",
x_dynlib)
system(cmd)
dyn.load(x_dynlib)
file.remove(c(x_mod, x_dynlib))
}
fortran_function <- function(..., signature, intent)
{
res <- .C64(...,
SIGNATURE = signature,
INTENT = intent,
NAOK = TRUE,
VERBOSE = 0)
return(res)
}
fortran_source("test.f90")
####################################################################################################
#import function
####################################################################################################
fortran_sum1d <- function(x)
{
res <- fortran_function("fortran_sum1d",
n = length(x), x = x, res = 0,
signature = c("int64", rep("double", 2)),
intent = c(rep("r", 2), "w"))$res
return(res)
}
fortran_sum2d <- function(x)
{
res <- fortran_function("fortran_sum2d",
n1 = nrow(x), n2 = ncol(x), x = x, res = 0,
signature = c(rep("int64", 2), rep("double", 2)),
intent = c(rep("r", 3), "w"))$res
return(res)
}
fortran_sum3d <- function(x)
{
res <- fortran_function("fortran_sum3d",
n1 = dim(x)[1], n2 = dim(x)[2], n3 = dim(x)[3], x = x, res = 0,
signature = c(rep("int64", 3), rep("double", 2)),
intent = c(rep("r", 4), "w"))$res
return(res)
}
fortran_lm <- function(x, y)
{
res <- fortran_function("fortran_lm",
n = min(length(x), length(y)), x = x, y = y, res = 0,
signature = c("int64", rep("double", 3)),
intent = c(rep("r", 3), "w"))$res
return(res)
}
fortran_nan_fill <- function(x, k)
{
res <- fortran_function("fortran_nan_fill",
n = length(x), x = x, k = k, res = numeric_dc(length(x)),
signature = c("int64", "double", "int64", "double"),
intent = c(rep("r", 3), "w"))$res
return(res)
}
####################################################################################################
#import data
####################################################################################################
set.seed(1)
test_vector1 <- rnorm(1e6)
test_vector2 <- rnorm(1e6)
test_matrix <- matrix(rnorm(1e6), nrow = 1e3)
test_array <- array(rnorm(1e6), dim = c(1e2, 1e2, 1e2))
####################################################################################################
#analysis
####################################################################################################
all.equal(sum(test_vector1), fortran_sum1d(test_vector1))
all.equal(sum(test_matrix), fortran_sum2d(test_matrix))
all.equal(sum(test_array), fortran_sum3d(test_array))
all.equal(lm(test_vector2 ~ test_vector1)$coefficients[[2]], fortran_lm(test_vector1, test_vector2))
all.equal(c(rep(NaN, 100), test_vector1[101:length(test_vector1)]),
fortran_nan_fill(test_vector1, 100))
####################################################################################################
Note: This tutorial presents an R-centered workflow emphasizing interactivity. Further customizations can be made to meet various needs and preferences. The R file needs to be run on RStudio and placed in the same directory as the Fortran file. Otherwise, the R and Fortran files’ path needs to be fully specified.