5

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?
4
  • 1
    Most often Fortran's swift Linear Algebra libraries are, in fact, BLAS (or one of its derivatives). Doesn't R use BLAS too ? Apr 15, 2015 at 9:29
  • A cursory search and it seems indeed that it does use those libraries, one still gets nowhere near the speed Fortran can achieve, as far as I can tell (anybody can dis/prove this?). I found this article benchmarking versions of R using different libraries, and it seems RRO using intel MKL is the winner. However, I would like this code to become a package eventually, so I have to stick with unmodified upstream R.
    – Ixxie
    Apr 15, 2015 at 10:19
  • The article you link to compares the performance of various implementations of BLAS. You would get very similar information from a comparison of Fortran linked to various implementations of BLAS. And the conclusion would be very similar -- on current architectures Intel's MKL is probably the fastest implementation of BLAS. Simply throwing Fortran in between R and slow BLAS won't change anything for the better. Apr 15, 2015 at 10:32
  • I should have been clearer: I am also aiming for portability. Changing the version of BLAS which R uses is thus not an option afaik; the acceleration should be done within a standard R installation. I am still not convinced Fortran will offer no advantages. Compilers may allow for optimizations that an interpreter cannot have, such as common subexpression elimination across multiple lines of code (common in ODE's). I do not know enough about this to be sure, but would like to test both options, so if anybody has some answer to my original question I am still interested in hearing it.
    – Ixxie
    Apr 15, 2015 at 11:52

2 Answers 2

3

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.

2
  • I've found Drew Schmidt's introductory packages very helpful too: github.com/wrathematics/Romp
    – Avraham
    Dec 4, 2018 at 2:12
  • There is the intrinsic iso_c_binding and bind() module in Fortran standard >2003 which can together give a Fortran subroutine/function the same feel of a C function. for example subroutine foo() bind(C,name="foo") can be called from C as void foo() without the need for the non-standard compiler-dependent trailing underscore "_".
    – Scientist
    Sep 13, 2019 at 19:47
3

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:

I am one of the authors of dotCall64 and spam.

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.