R: Advantages of using a Fortran subroutine with .Call and C/C++ wrapper instead of .Fortran?

Question

I have an R package which uses lots of Fortran subroutines for nested loops of recursive linear algebra computations (depending heavily on BLAS and LAPACK routines). As an interface to Fortran, I use .Fortran function. I just read Jonathan Callahan's blog post about using .Call instead of .C in case of subroutines written in C/C++, and it got me thinking that would it be better to use .Call interface also when using Fortran subroutines, by writing a simple wrapper in C which then calls the Fortran subroutines?

As said, my Fortran codes are quite simple in a sense that I just play with multidimensional arrays of type double or integer. But I have learned that I must write quite a lot of checks in R side to ensure that everything doesn't crash because of I accidentally forgot to change the storage mode of some matrix to integer or the dimensions of some matrix were changed etc.

Subroutines are written as F90/95.

Sloop

seems reasonable to use .Call() with some C function, and you can then indeed call your Fortran subroutines from the C code, which is relatively easy (or even do everything in C if you do not really require Fortran). – Whitehead

Yes but what kind of benefits would that bring if any? I could switch to C completely, but I that would be too cumbersome and I doubt it would be useful as I would then be calling Fortran BLAS functions from C anyway. – Sloop

Might be relevant: ualberta.ca/AICT/RESEARCH/LinuxClusters/doc/ifc91/main_for/… – Someway

Why would you switch here? You are then will have another layer of abstraction. Calling a C routine which in turn calls FORTRAN can never be quicker than what you have, R -> FORTRAN. – Seventh

I don't think it's that simple, my understanding based on the link on my question and other sources, is that the overheads of .C and .Fortran are much larger than .Call at least because the extra copying of the arguments. There might be some performance differences with the type checking in C compared to R etc also, don't know about that. – Sloop

Answer 1

There might be an advantage if you are working with a large dataset. .Call can be much faster because you are not copying the data each time you call the function. For the case described in this question, there will be no such advantage, because the R 2.15.1 release notes state

.C() and .Fortran() do less copying: arguments which are raw, logical, integer, real or complex vectors and are unnamed are not copied before the call, and (named or not) are not copied after the call. Lists are no longer copied (they are supposed to be used read-only in the C code).

Switching to .Call means you give up the convenience of the .Fortran interface. You'd pass SEXPs into the C code, do any checks/manipulation of the data using the (scary and not well-documented) R API, and then call a Fortran function from C. Anyone else working with your code will have to understand the R API and C/Fortran interop.

Aemia

Yes the copying seems to be one clear benefit, although as of R 2.15.1 the difference is slightly reduced: ".C() and .Fortran() do less copying: arguments which are raw, logical, integer, real or complex vectors and are unnamed are not copied before the call, and (named or not) are not copied after the call." (From R News). – Sloop

It looks like in my application the effects of copying are not significant; I tested calling my code with .Fortran using DUP=FALSE and DUP=TRUE, and saw practically no differences at all, although I was using arrays containing millions of elements. – Sloop

I was not aware of the change in behavior of those functions. I've edited the question to reflect it. – Aemia

Answer 2

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.

Vola