Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Efficient calling of F95 in R: use .Fortran or .Call?

Tags:

c

r

dll

fortran

I am writing some R simulation code, but want to leverage Fortran's swift Linear Algebra libraries to replace the core iteration loops. So far I was looking primarily at the obvious option of using .Fortran to call linked F95 subroutines; I figured I should optimize memory use (I am passing very large arrays) and set DUP=FALSE but then I read the warning in manual about the dangers of this approach and of its depreciation in R 3.1.0 and disablement in R 3.2.0. Now the manual recommends switching to .Call, but this function offers no Fortran support natively.

My googling has yielded a stackoverflow question which explores an approach of linking the Fortran subroutine through C code and the calling it using .Call. This seems to me the kind of thing that could either work like a charm or a curse. Hence, my questions:

  1. Aiming for speed and robustness, what are the risks and benefits of calling Fortran via .Fortran and via .Call?
  2. Is there a more elegant/efficient way of using .Call to call Fortran subroutines?
  3. Is there another option altogether?
like image 279
Ixxie Avatar asked Apr 15 '15 09:04

Ixxie


2 Answers

Here's my thoughts on the situation:

.Call is the generally preferred interface. It provides you a pointer to the underlying R data object (a SEXP) directly, and so all of the memory management is then your decision to make. You can respect the NAMED field and duplicate the data if you want, or ignore it (if you know that you won't be modifying the data in place, or feel comfortable doing that for some other reason)

.Fortran tries to automagically provide the appropriate data types from an R SEXP object to a Fortran subroutine; however, its use is generally discouraged (for reasons not entirely clear to me, to be honest)

You should have some luck calling compiled Fortran code from C / C++ routines. Given a Fortran subroutine called fortran_subroutine, you should be able to provide a forward declaration in your C / C++ code as e.g. (note: you'll need a leading extern "C" for C++ code):

void fortran_subroutine_(<args>);

Note the trailing underscore on the function name -- this is just how Fortran compilers (that I am familiar with, e.g. gfortran) 'mangle' symbol names by default, and so the symbol that's made available will have that trailing underscore.

In addition, you'll need to make sure the <args> you choose map to from the corresponding C types to the corresponding Fortran types. Fortunately, R-exts provides such a table.

In the end, R's R CMD build would automatically facilitate the compilation + linking process for an R package. Because I am evidently a glutton for punishment, I've produced an example package which should provide enough information for you to get a sense of how the bindings work there.

like image 165
Kevin Ushey Avatar answered Oct 17 '22 00:10

Kevin Ushey


3. Is there another option altogether? YES

The R package dotCall64 could be an interesting alternative. It provides .C64() which is an enhanced version of the Foreign Function Interface, i.e., .C() and .Fortran().

The interface .C64() can be used to interface both Fortran and C/C++ code. It

  • has a similar usage as .C() and .Fortran()
  • provides a mechanism to avoid unnecessary copies of read-only and write-only arguments
  • supports long vectors (vectors with more then 2^31-1 elements)
  • supports 64-bit integer type arguments

Hence, one can avoid unnecessary copies of read-only arguments while avoiding the .Call() interface in combination with a C wrapper function.

Some links:

  • The R package dotCall64: https://CRAN.R-project.org/package=dotCall64
  • Description of dotCall64 with examples: https://doi.org/10.1016/j.softx.2018.06.002
  • An illustration where dotCall64 is used to make the sparse matrix algebra R package spam compatible with huge sparse matrices (more than 2^31-1 non-zero elements): https://doi.org/10.1016/j.cageo.2016.11.015

I am one of the authors of dotCall64 and spam.

like image 26
Nairolf Avatar answered Oct 17 '22 00:10

Nairolf