Passing multidimensional array from R to Fortran using Rcpp

Question

I'm new to C++ and the usage of Rcpp in R. I have many Fortran subroutines that I would like to include in an R package. I have followed this example, and it works perfectly. Now, I would like to implement one of my subroutines, a dummy Fortran subroutine saved in the file mytest.f90, and it looks like this:

module mytest
    use, intrinsic :: iso_c_binding, only: c_int

    implicit none

    public :: mytest_f

contains

    subroutine mytest_f(l, m, n, a, c) bind(C)

    integer(kind=c_int) :: l, m, n
    integer(kind=c_int) :: a(l, m, n), c(l,m,n)
    c = c - a
    end subroutine mytest_f
end module mytest

The Rcpp code and the C-wrapper are saved in one file called Rcppmytest.cpp and it contains:

#include <Rcpp.h>
using namespace Rcpp;

//C wrapper ??
extern "C" {
  void mytest_f(int l, int m, int n, int a, int c);
}

// Rcpp
int c_mytest_f(Rcpp::IntegerVector a3D) {

  // receives the 3D array
  Rcpp::IntegerVector DIM = a3D.attr("dim");

  int l, m, n;
  l = DIM[0];
  m = DIM[1];
  n = DIM[2];

  int c[l][m][n];

  mytest_f(l, m, n, a3D, c);

  return c;
}

However, I have two errors in the IDE (i.e. Rstudio). Firstly, in line mytest_f(l, m, n, a, c);, and the error is

not matching function for call to "mytest_f"

While the second error is:

Error: Cannot initialize return object of type "int *" with an lvalue of type 'int [a_size]'

When I try to get more info from the compilation process using devtools::document() in R, this error message appears:

   Rcppmytest.cpp:n function ‘int c_mytest_f(Rcpp::IntegerVector)
   Rcppmytest.cpp:21:7:warning: O C++ forbids variable length array ‘c-Wvla
      21 |   int cze];
         |       ^
   Rcppmytest.cpp:22:7:warning: O C++ forbids variable length array ‘a-Wvla
      22 |   int aze];
         |       ^
   Rcppmytest.cpp:35:21:error: valid conversion from ‘int*to ‘int-fpermissive
      35 |   mytest_f(l, m, n, a
         |                     ^
         |                     |
         |                     int*
   Rcppmytest.cpp:6:42:note: initializing argument 4 of ‘void mytest_f(int, int, int, int, int)
       6 |   void mytest_f(int l, int m, int n, int a c);
         |                                      ~~~~^
   Rcppmytest.cpp:35:24:error: valid conversion from ‘int*to ‘int-fpermissive
      35 |   mytest_f(l, m, n, a, c
         |                        ^
         |                        |
         |                        int*
   Rcppmytest.cpp:6:49:note: initializing argument 5 of ‘void mytest_f(int, int, int, int, int)
       6 |   void mytest_f(int l, int m, int n, int a, int c
         |                                             ~~~~^
   Rcppmytest.cpp:37:10:error: valid conversion from ‘int*to ‘int-fpermissive
      37 |   return c
         |          ^
         |          |
         |          int*
   make: *** [/usr/lib/R/etc/Makeconf:178: Rcppmytest.o] Error 1
   ERROR: compilation failed for package ‘mytestpackage’

My questions would be:

1. How to make it work (Any help and hints would be appreciated)?

2. Should I care about the column-/row-major order difference between R/Fortran and C/C++?

Notes: The reason I'm using Rcpp is that I tried writing a small init.c file similar to here however, my efforts wasn't successful.

Answer 1

You need to adjust the C++ side of your glue, and the good news is that it is simple -- at least to get to a pointer and pass it along. Here is what I would do:

#include <Rcpp.h>

extern "C" {
    void mytest_f(int l, int m, int n, int *a, int *c);
}

// [[Rcpp::export]]
Rcpp::IntegerVector c_mytest_f(Rcpp::IntegerVector a3D) {

    // receives the 3D array
    Rcpp::IntegerVector dim = a3D.attr("dim");

    int l = dim[0],
        m = dim[1],
        n = dim[2];

    Rcpp::IntegerVector c(l*m*n); // vectors are contiguous

    mytest_f(l, m, n, &(a3D[0]), &(c[0]));

    return c;
}

The changes are to

• pass vectors a and c as pointers to double along with the dimension information
• allocate the result vector as another Rcpp::IntegerVector which gets you both the allocation and the return object
• pass a and c as pointers by accessing the address of their first element.

Lastly, files like the init.c you mention are autowritten for you (at least for the C / C++ part) and as mytest_f is not called from R you need to worry about it in there -- the calling is done by the glue code, and Rcpp makes setting up glue fairly easy.

All that said, you may have a problem with row-versus-column order between C (or C++) and Fortran. In the 2d-case you could do what you did here with RcppArmadillo types and transpose, I am not aware of a 'quick' fix for the 3d case though there may be a common trick.