Yutaka Masuda has a great tutorial of calling Fortran subroutines from R with its Foreign Function Interface (FFI).
https://masuday.github.io/fortran_tutorial/r.html
This FFI approach has two limitations. First, it copies R objects, which hurts performance when passing a large R object to Fortran subroutines. 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 very long vector or large matrix/array to Fortran subroutines because its dimension needs to be an integer. Fortunately, We can use the R package dotCall64 to deal with these limitations.
Below is a Fortran file named example.f containing subroutines to be called from R. These subroutines are used to illustrate how to pass a vector, a matrix, an array (rank higher than two), and multiple vectors, from R to Fortran. OpenMP is used to showcase Fortran’s parallel computing feature.
!example.f
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine vector_sum(n, x, res)
integer(8) :: n, i
real(8) :: x(n), res
res = 0d0
!$omp parallel do reduction(+:res)
do i = 1, n
res = res + x(i)
end do
end subroutine
subroutine matrix_sum(n, m, x, res)
integer(8) :: n, m, i, j
real(8) :: x(n, m), res
res = 0d0
!$omp parallel do reduction(+:res) collapse(2)
do j = 1, m
do i = 1, n
res = res + x(i, j)
end do
end do
end subroutine
subroutine array_sum(n1, n2, n3, x, res)
integer(8) :: n1, n2, n3, i, j, k
real(8) :: x(n1, n2, n3), res
res = 0d0
!$omp parallel do reduction(+:res) collapse(3)
do k = 1, n3
do j = 1, n2
do i = 1, n1
res = res + x(i, j, k)
end do
end do
end do
end subroutine
subroutine vector_lm(n, x, y, res)
integer(8) :: n, i
real(8) :: x(n), y(n), res
real(8) :: sum_x, sum_y, sum_xx, sum_xy
sum_x = 0d0
sum_y = 0d0
sum_xx = 0d0
sum_xy = 0d0
!$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 subroutine
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Below is an R file calling the above subroutines and validating the results with corresponding R codes.
#header
####################################################################################################
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
fortran_source <- function(x)
{
#x: fortran file
x_dynlib <- paste0(substr(x, 1, nchar(x) - 2), .Platform$dynlib.ext)
cmd <- paste("gfortran",
"-ffree-form",
"-fimplicit-none",
"-fopenmp",
"-fpic",
"-shared",
x,
"-o",
x_dynlib)
system(cmd)
dyn.load(x_dynlib)
file.remove(x_dynlib)
}
fortran_function <- function(..., signature, intent)
{
res <- dotCall64::.C64(...,
SIGNATURE = signature,
INTENT = intent,
NAOK = TRUE,
VERBOSE = 0)
return(res)
}
fortran_source("example.f")
####################################################################################################
#import function
####################################################################################################
fortran_vector_sum <- function(x)
{
res <- fortran_function("vector_sum",
n = length(x), x = x, res = 0,
signature = c("int64", rep("double", 2)),
intent = c(rep("r", 2), "w"))$res
return(res)
}
fortran_matrix_sum <- function(x)
{
res <- fortran_function("matrix_sum",
n = nrow(x), m = ncol(x), x = x, res = 0,
signature = c(rep("int64", 2), rep("double", 2)),
intent = c(rep("r", 3), "w"))$res
return(res)
}
fortran_array_sum <- function(x)
{
res <- fortran_function("array_sum",
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_vector_lm <- function(x, y)
{
res <- fortran_function("vector_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)
}
####################################################################################################
#import data
####################################################################################################
set.seed(1)
test_vector1 <- rnorm(1e6)
test_vector2 <- rnorm(1e6)
test_matrix <- matrix(rnorm(1e6),
nrow = 1e2)
test_array <- array(rnorm(1e6), dim = c(1e2, 1e2, 1e2))
####################################################################################################
#analysis
####################################################################################################
all.equal(sum(test_vector1), fortran_vector_sum(test_vector1))
all.equal(sum(test_matrix), fortran_matrix_sum(test_matrix))
all.equal(sum(test_array), fortran_array_sum(test_array))
all.equal(lm(test_vector2 ~ test_vector1)$coefficients[[2]], fortran_vector_lm(test_vector1, test_vector2))
####################################################################################################
Note: This tutorial presents an R-centered workflow emphasizing interactivity. Further customizations can be made to meet various demands and preferences. To ensure that the Fortran subroutines work with a very long vector or large matrix/array passed from R, 64-bit integers are used. 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.