Call parallel fortran MPI subroutine from R

Question

I would like to write some parallel Fortran code in a subroutine that can be called by R (I would like to read in data from R and send it to a parallel Fortran MPI). I have noticed, however, that when I run the following program as a subroutine (i.e. substitute "subroutine" for "program"), the code no longer compiles (it does compile when it is a program). I am compiling the code using the mpif90 from MPICH in Linux.

Is it possible to initialize and finalize an MPI in a subroutine in Fortran? If not, is it still possible to somehow call a parallel Fortran MPI from R? If not in Fortran, can this be done in C?

Here's the code:

module global
  integer numnodes,myid,mpi_err
  integer, parameter :: my_root=0
end module global

module fmpi
  include 'mpif.h'
end module fmpi

subroutine init
  use fmpi
  use global
  implicit none
  call MPI_INIT( mpi_err )
  call MPI_COMM_SIZE( MPI_COMM_WORLD, numnodes, mpi_err )
  call MPI_Comm_rank(MPI_COMM_WORLD, myid, mpi_err)
end subroutine init

program test
  use global
  use fmpi
  implicit none
  real*8:: dat(10)
  integer*4:: i
  call init
  if(myid == my_root) then
    do i=1,10
      dat(i) = i
    enddo
    print *,dat(1)
  endif
  call mpi_finalize(mpi_err)
end program test

Answer 1

Here is a simple Fortran/MPI subroutine that I want to call from R:

subroutine test(id, ierr)
  use mpi
  implicit none
  integer*4 id, ierr
  call MPI_Comm_rank(MPI_COMM_WORLD, id, ierr)
end subroutine test

To call this from R on a Linux machine, I built a shared object file using the Open MPI wrapper command "mpif90":

$ mpif90 -fpic -shared -o test.so test.f90

I tried to use "R CMD SHLIB", but eventually decided that it was easier to get "mpif90" to create a shared object than to get "R CMD SHLIB" to deal with MPI. The downside is that the command is gfortran specific. For a different compiler, you might get some help by using the "SHLIB" --dry-run option:

$ R CMD SHLIB --dry-run test.f90

This will display the commands that it would have used to create the shared object using your compiler. You can then modify the commands to use "mpif90" in order to handle the MPI headers and libraries.

Here is an R script that calls the Fortran test subroutine. It loads Rmpi (which automatically calls MPI_Init), loads the shared object containing my Fortran subroutine, and then calls it:

# SPMD-style program: start all workers via mpirun
library(Rmpi)
dyn.load("test.so")

# This Fortran subroutine will use MPI functions
r <- .Fortran("test", as.integer(0), as.integer(0))

# Each worker displays the results
id <- r[[1]]
ierr <- r[[2]]
if (ierr == 0) {
  cat(sprintf("worker %d: hello\n", id))
} else {
  cat(sprintf("ierr = %d\n", ierr))
}

# Finalize MPI and quit
mpi.quit()

Since it's an SPMD-style program, it doesn't spawn workers, like many Rmpi examples. Instead, all of the workers are started via mpirun, which is the typical way of executing C and Fortran MPI programs:

$ mpirun -n 3 R --slave -f test.R

This runs three instances of my R script, so the output is:

worker 0: hello
worker 1: hello
worker 2: hello

I think that structuring the code in this way makes it easy to use MPI from R and from any number of Fortran subroutines.