Search code examples
rcmdfortran

Using blas/Lapack in Fortran file for R


I have the following problem. I have a file say file.f90 in which I have implemented some Fortran subroutines, say called foo. Then I compile these functions using "R CMD SHLIB file.f90". To use the function foo in a seperate R file I then use dyn.load("foo.dll") and to call it .Fortran("foo", ...).

So far so good. But now I need to use some functions implemented in Lapack.

I have no idea how to do this or where to have a look. I have only tried calling "R CMD SHLIB file.f90 -llapack" but already there I get an error that llapack has not been found. Any hints would be greatly appreciated!!

EDIT: I have finally found an answer to my question with the help of everyone here and with looking up much on the internet. I have to say the solution is quite easy but as I am quite a noob when it comes to these things it still took some time. So here s my solution for Windows 11 and R studio 4.1

Assume that our R session/project has the path PATH_PROJ, e.g "C:\Users\Myname\Documents\MyProject". Then I created a new folder named "f90files" in which I intended to save all Fortran functions, so PATH_PROJ\f90files. Next, I needed the path of my R's Lapack PATH_LAPACK, e.g "C:\Program Files\R\R-4.1.2\bin\x64\Rlapack.dll".

In PATH_PROJ\f90files I then implemented the Fortran subroutine as suggested by Jean-Claude Arbaut:

subroutine eigvals(n, a, vre, vim, info)
    implicit none
    integer :: n, info
    integer, parameter :: lwork = 65536
    double precision :: a(n, n), vre(n), vim(n)
    double precision, save :: work(lwork)

    call dgeev("n", "n", n, a, n, vre, vim, 0d0, 1, 0d0, 1, work, lwork, info)
end subroutine

Following this, I started up Windows command prompt and typed

gfortran -shared PATH_LAPACK PATH_PROJ\f90files\eigvals.f90 -o PATH_PROJ\f90files\eigvals.so -o PATH_PROJ\f90files\eigvals.dll

and further

gfortran -shared PATH_LAPACK PATH_PROJ\f90files\eigvals.f90 -o PATH_PROJ\f90files\eigvals.so

(maybe this can be done in one go?)

With this all was nicely compiled. In R I then loaded the function using

dyn.load("PATH_PROJ\f90files\eigvals.dll")

Finally, using the implementation given below, I ran

eigvals <- function(a) {
  if (is.matrix(a) && is.double(a) && nrow(a) == ncol(a)) {
    n <- nrow(a)
    s <- .Fortran("eigvals", n = as.integer(n), a = a, vre = double(n), vim = double(n), info = 0L)
    structure(complex(real = s$vre, imaginary = s$vim), info = s$info)
  } else stop("Invalid input")
}


eigvals(a)

and voilà we are done! Thanks again to everyone!


Solution

  • The libraries you are looking for are Rblas.dll and Rlapack.dll in the R-4.2.2\bin\x64 directory (replace 4.2.2 with your version).

    Here is an example. Let's compute eigenvalues using LAPACK's dgeev.

    Fortran file eigvals.f90. Here to simplify lwork is a constant, but in "real" code you would have to do this more carefully.

    subroutine eigvals(n, a, vre, vim, info)
        implicit none
        integer :: n, info
        integer, parameter :: lwork = 65536
        double precision :: a(n, n), vre(n), vim(n)
        double precision, save :: work(lwork)
    
        call dgeev("n", "n", n, a, n, vre, vim, 0d0, 1, 0d0, 1, work, lwork, info)
    end subroutine
    

    Compile with either one of the following commands (change the path as necessary). If you are on Windows, do this from the Rtools bash window. On Linux the extension is .so and not .dll.

    gfortran -shared -L/c/App/R/R-4.2.2/bin/x64 eigvals.f90 -lRlapack -o eigvals.dll
    R CMD SHLIB eigvals.f90 -lRlapack
    

    In R, you can now do, assuming you have setwd() to the directory containing the DLL:

    a <- matrix(c(2, 9, 4, 7, 5, 3, 6, 1, 8), 3, 3, byrow = T)
    
    dyn.load("eigvals.dll")
    is.loaded("eigvals")
    
    eigvals <- function(a) {
      if (is.matrix(a) && is.double(a) && nrow(a) == ncol(a)) {
        n <- nrow(a)
        s <- .Fortran("eigvals", n = as.integer(n), a = a, vre = double(n), vim = double(n), info = 0L)
        structure(complex(real = s$vre, imaginary = s$vim), info = s$info)
      } else stop("Invalid input")
    }
    
    eigvals(a)
    

    Pay attention to the number and type of subroutine arguments in .Fortran, otherwise you may crash R. Note also that you must call Fortran subroutines, not functions.

    I you are on Windows and using R-4.2 with Rtools 4.2, there is an extra trick: the compiler is no longer in the default directories. See this. You have first to do, in the Rtools bash window:

    export PATH=/x86_64-w64-mingw32.static.posix/bin:$PATH
    

    If you are using the compiler from the Windows command prompt, you will have to modify the PATH environment variable accordingly.