Fortran calling R

1

Does anyone call R from Fortran?

R and its packages have almost every statistics algorithm I want, and I would like to access it from Fortran. A simple but crude way is to write unformatted stream data from Fortran, which R reads and processes. For example, the code

program xcall_r
! 12/25/2020 10:12 AM demonstrate Fortran calling R, passing data with binary file
implicit none
integer :: iter
integer, parameter :: n = 1000000,bin_unit = 20, niter = 10
real(kind=kind(1.0d0)) :: x(n)
character (len=*), parameter :: bin_file = "double.bin"
logical, parameter :: call_r = .true.
do iter=1,niter
   call random_number(x)
   write (*,"(a,f11.7)") "mean = ",sum(x)/n
   if (call_r) then
      open (unit=bin_unit,file=bin_file,action="write",access="stream",form="unformatted",status="replace")
      write (bin_unit) x
      close (bin_unit)
      call execute_command_line("C:\programs\R\R-4.0.1\bin\x64\rterm.exe --vanilla --slave  < xread_bin.r")
   end if
end do
end program xcall_r

for R script xread_bin.r

inp = file("double.bin","rb")
x = readBin(inp, "double",n=10000000) # n is max number of values to read -- can read fewer
cat("from R: ",mean(x),"\n\n")

gives output such as

mean =   0.5000192
from R:  0.5000192 

mean =   0.5003266
from R:  0.5003266 
<snip>

The wall time elapsed is 2.2 s when the program calls R 10 times and 0.1 s when it does not, so calling R has overhead of about 0.2 s (even if R script does nothing). Passing estimated parameters from R back to Fortran instead of just printing them would be more involved. It is possible to call C++ functions from R using the Rcpp package. There is an RFortran library, but it works only with Intel Fortran on Windows.

2

If you look at section 5.11 there is the possibility to evaluate R expressions in C. But the R API is based quite heavily on macros and define statements. It will likely take a lot of work to get something intuitive and simple to use.

3

Apparently, you can link directly to the R runtime and create an embedded R environemnt. But you will need to write some C code.

4

FWIW, I once had to compile R sources and noticed that the compilation required a Fortran compiler. Turns out that R uses Fortran implementations of various statistical routines. I don’t know the exact routines you need, but it may just be possible that the functionality you need is implemented in Fortran in the first place. In this case, rather than Fortran -> R -> Fortran, it may be possible to remove R out of the picture. :grinning:

2 Likes
5

Sounds it bit like “calling Excel from R”
As @eugene_epshteyn pointed out many R packages are in C or FORTARN, and personally I have embedded some of them directly in my libraries. However, many of the quite advanced stat packages have a large body of R code and the contained C and Fortran routines are just for the heavy lifting but not to postprocess the results as you finally see them in R. Maybe you have a lock at the intel mkl vsl and intel oneDal.

7

I love R and use it daily. But normally I call compiled languages from R, not the other way around. However, there is RInside. It’s by the author of Rcpp, Dirk Eddelbuettel

2 Likes
8

@Beliavsky, @edsterjo,

I have created a small example using RInside to call R from Fortran at: GitHub - ivan-pi/Fortran-RInside: Demonstration of Fortran interface to Rinside

Adapting your example it becomes as simple as:

program xcall_r
! Demonstrate Fortran calling R, passing data with binary file
use RInside_interface
implicit none
integer :: iter
integer, parameter :: n = 100000, bin_unit = 20, niter = 3
real(kind=kind(1.0d0)) :: x(n)
character (len=*), parameter :: bin_file = "double.bin"
logical, parameter :: call_r = .true.

call setupRinC()

do iter = 1, niter
  call random_number(x)
  write (*,"(a,f11.7)") "mean = ",sum(x)/n
  if (call_r) then
    open (unit=bin_unit,file=bin_file,action="write",access="stream",form="unformatted",status="replace")
    write (bin_unit) x
    close (bin_unit)
    call xread_bin()
  end if
end do

call teardownRinC()

contains

  subroutine xread_bin()
    type(SEXP) :: res
    type(SEXP) :: inp, x
    res = evalInR('inp = file("double.bin","rb")')

    res = evalInR('x = readBin(inp, "double",n=100000)') ! n is max number of values to read -- can read fewer
    res = evalInR('cat("from R: ",mean(x),"\n\n")')
    res = evalInR('close(inp)')
  end subroutine

end program xcall_r

and produces the output:

mean =   0.4990118
from R:  0.4990118 

mean =   0.4982012
from R:  0.4982012 

mean =   0.4998608
from R:  0.4998608

For some reason if the array is too large (1000000 values) I get errors such as:

Error: C stack usage  8009524 is too close to the limit
Error: C stack usage  8009572 is too close to the limit
Error: C stack usage  8009460 is too close to the limit
Fatal error: unable to initialize the JIT

which seems to be related to the stack size settings of the R binary.

2 Likes
9

Thanks, it works for me on Windows Subsystem for Linux. (I wonder how to modify the make file to work on Windows). A convenient alternative to running evalInR on each line of R code is to run evalInR on (source(“script.r”)), where script.r contains the R commands, although there is a 50% speed penalty in my tests. The slightly modified Fortran code is

program xcall_r
! Demonstrate Fortran calling R, passing data with binary file
use RInside_interface, only: evalInR, SEXP, teardownRinC, setupRinC
implicit none
integer :: iter,t1,t2,itick
integer, parameter :: n = 100000, bin_unit = 20, niter = 100
real(kind=kind(1.0d0)) :: x(n)
character (len=*), parameter :: bin_file = "double.bin"
logical, parameter :: call_r = .true.
call system_clock(t1)
call setupRinC()

do iter = 1, niter
  call random_number(x)
  write (*,"(a,f11.7)") "mean = ",sum(x)/n
  if (call_r) then
    open (unit=bin_unit,file=bin_file,action="write",access="stream",form="unformatted",status="replace")
    write (bin_unit) x
    close (bin_unit)
    call xread_bin()
  end if
end do

call teardownRinC()
call system_clock(t2,itick)
write (*,"(/,a,f0.4)") "time elapsed (s) = ",(t2-t1)/dble(itick)
contains

  subroutine xread_bin()
    type(SEXP) :: res
    type(SEXP) :: inp, x
    logical, parameter :: call_script = .false.
    if (call_script) then
       res = evalinR('source("print_mean.r")')
    else
       res = evalInR('inp = file("double.bin","rb")')
       res = evalInR('x = readBin(inp, "double",n=100000)') ! n is max number of values to read -- can read fewer
       res = evalInR('cat("from R: ",mean(x),"\n\n")')
       res = evalInR('close(inp)')
    end if
  end subroutine xread_bin

end program xcall_r

which calls print_mean.r

inp = file("double.bin","rb")
x = readBin(inp,"double",n=100000)
cat("from R: ",mean(x),"\n\n")
close(inp)
10

Unfortunately, I have very little experience using Fortran in Windows. Some searching suggests that as long as you are using the cygwin environment the Makefile should work per usual.

I wrapped additionally the function Rf_allocVector to create a R vector of doubles directly in Fortran. The function passToR is then used to pass it directly to the R session:

type(SEXP), target :: x_R
real(kind=kind(1.0d0)), pointer :: x_ptr(:) => null()
integer, parameter :: REALSXP = 14
type(SEXP) :: dummy

call setupRinC()

! Allocate R vector of reals object
x_R = Rf_protect(Rf_allocVector(REALSXP,n))

! Recover pointer to the underlying object memory
x_ptr => double_from_SEXP(x_R,n)

do iter = 1, niter
  call random_number(x_ptr)

  write (*,"(a,f11.7)") "mean = ",sum(x_ptr)/n
  
  call passToR(x_R,'x')
  dummy = evalInR('cat("from R: ",mean(x),"\n\n")')

end do

call Rf_unprotect(1)

call teardownRinC()

Edit: the code is in the branch GitHub - ivan-pi/Fortran-RInside at alloc.

It should be possible to retrieve the mean value calculated in R, by assigning the result of evalInR("mean(x)") to a new SEXP, and using one of the macros in Rinternals.h to recover it’s value as a scalar double.

Edit2: to recover the double precision value the function Rf_asReal(x) is available. Instead of printing the mean in R, you can assign the value to a SEXP and retrieve the value as follows:

  call passToR(x_R,'x')
  res = evalInR("mean(x)")
  write(*,"(a,f11.7)") "from R: ", Rf_asReal(res)

It might be nicer to wrap the C++ interface of RInside, but this would involve more work. I have the feeling the C++ interface does not require the protect and unprotect calls, and manages these through the class constructors and destructors.