Extend R with Fortran

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.

6 Likes

Seems good, bookmarked.

The R package sparsegl uses dotCall64 to call Fortran from R.

The most popular R IDE RStudio is scheduled to add Fortran syntax highlighting in 2023.

2 Likes