Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

R .C interface: Passing multidimensional arrays

I've written up a function "foo" in C that I want to call from an R program. The function takes in as input a matrix and does some operations on it, (say add 1 to each element). While it is easy to a single vector as

.C("foo", n=as.integer(5), x=as.double(rnorm(5)))

with foo implemented as

void foo(int *nin, double *x)
{
int n = nin[0];

int i;

for (i=0; i<n; i++)
    x[i] = x[i] * x[i];
} 

How do I pass in a two dimensional array? If I change the "double *x" to "double **x" it gives a segmentation fault. Any pointers appreciated.

like image 829
broccoli Avatar asked Oct 22 '12 02:10

broccoli


2 Answers

No need to give up on .C for straight-forward manipulations like this. Remember that a matrix in R is a vector + dimensions. Likewise in C, so pass the matrix and its dimensions, and access elements of the matrix as appropriate offsets into a vector. Something like

void cplus1(double *x, int *dim)
{
    for (int j = 0; j < dim[1]; ++j)
        for (int i = 0; i < dim[0]; ++i)
            x[j * dim[0] + i] += 1;
}

so using inline as a nice party trick

library(inline)
sig <- signature(x="numeric", dim="integer")
body <- "
    for (int j = 0; j < dim[1]; ++j)
        for (int i = 0; i < dim[0]; ++i)
            x[j * dim[0] + i] += 1;
"

cfun <- cfunction(sig, body=body, language="C", convention=".C")

plus1 <- function(m) {
    m[] = cfun(as.numeric(m), dim(m))$x
    m
}
like image 190
Martin Morgan Avatar answered Oct 20 '22 04:10

Martin Morgan


Give up on .C() and switch to .Call() which lets you pass whole R objects as so-called SEXP objects.

You can parse those the hard way via the C API of R, or via Rcpp with (what we think) nice higher-level abstractions.

R> library(inline)  # use version 0.3.10 for rcpp() wrapper
R> 
R> addMat <- rcpp(signature(ms="numeric"), body='
+    Rcpp::NumericMatrix M(ms);
+    Rcpp::NumericMatrix N = Rcpp::clone(M);
+    for (int i=0; i<M.nrow(); i++)
+       for (int j=0; j<M.ncol(); j++) 
+          N(i,j) = M(i,j) + 1;
+    return(N);
+ ')
R> addMat(matrix(1:9,3,3))
     [,1] [,2] [,3]
[1,]    2    5    8
[2,]    3    6    9
[3,]    4    7   10
R> 
like image 20
Dirk Eddelbuettel Avatar answered Oct 20 '22 02:10

Dirk Eddelbuettel