Introduction to Calling C and Fortran from R

1. Introduction to Calling C and Fortran from R In UNIX

2. Why Call C or Fortran from R? • R is rotten at iterative algorithms that require loops iterated many times (especially Markov chain Monte Carlo), not as bad as S-PLUS, but still bad. A way to get all the speed advantages of C or Fortran with most of the convenience of R is to write the inner loop in C and call it from R.

3. C functions and Fortran subroutines callable from R C functions: void foo(long *nin, double *x) { int n = nin[0]; int i; for (i=0; i<n; i++) x[i] = x[i] * x[i]; }

4. two properties that are required for a function callable from R • The function does not return a value. All work is accomplished as a "side effect" (changing the values of arguments). • All the arguments are pointers. Even scalars are vectors (of length one) in R.

5. Fortran subroutines subroutine bar(n, x) integer n double precision x(n) integer I do 100 i = 1, n x(i) = x(i) ** 2 100 continue end

6. Compiling and Dynamic Loading • compile it to a shared library. The command R CMD SHLIB foo.c • Now the code can be dynamically loaded into R by doing (in R) dyn.load("foo.so")

7. The Call to C • The actual call to C from R is made by the R function .C like this .C("foo", n=as.integer(5),x=as.double(rnorm(5)))

8. The function .C takes nargs + 1 arguments if the C function being called takes nargs arguments. • The first argument is the name of the C function to be called. • the rest are the arguments to the function call. The arguments generally have to be coerced to the proper type using as.integer, as.double, and similar functions. The function .C returns a list whose elements are the arguments to the function call (in this case n and x). The values are the values after the call (as changed by the function foo).

9. The Fortran is similar. In UNIX do R CMD SHLIB bar.f and in R do dyn.load("bar.so") .Fortran("bar", n=as.integer(5), x=as.double(rnorm(5))) We just change the call from .C to .Fortran

10. R Wrapper Functions foo <- function(x) { if (!is.numeric(x)) stop("argument x must be numeric") out <- .C("foo", n=as.integer(length(x)), x=as.double(x)) return(out$x) }

11. R Wrapper Functions has three benefits • It allows some error checking in R, where it is easier than in C. • It allows some arguments (like n here) to be calculated so they don't have to be supplied by the user. • It allows you to return only what the user needs. Now the call is much simpler, just foo(x)

12. Using R Random Number Generators • The R random number generators and also all the other functions for probability distributions (not only rnorm but also dnorm, pnorm, and qnorm and so forth for other distributions) are callable from C (but not Fortran).

13. example #include <R.h> #include <Rmath.h> void baz(int *nin, double *sin, double *tin, double *x) { int n = nin[0]; double s = sin[0]; double t = tin[0]; int i; GetRNGstate(); for (i = 0; i < n; i++) x[i] = rbeta(s, t); PutRNGstate(); } GetRNGstate(); PutRNGstate(); must be done before and after (respectively) all calls to random number functions (in this example rbeta).

14. set.seed(42) rbeta(5, 1.5, 2.5) dyn.load("baz.so") baz <- function(n, s, t) { .C("baz", n = as.integer(n), s = as.double(s), t = as.double(t), x = double(n))$x } set.seed(42) baz(5, 1.5, 2.5)

15. Error Messages from C Use the R error or warning function. They work just like printf but produce an R error or warning as the case may be. Here is an example. #include <R.h> void qux(int *nin, double *x) { int n = nin[0]; int i; if (n < 1) error("arg n was %d, must be positive\n", n); for (i = 0; i < n; i++) x[i] = x[i] * x[i]; } And try it out. dyn.load("qux.so") .C("qux", n = as.integer(0), x = as.double(1:4))

16. Motivation: • Speed • Efficient memory management • Using existing C libraries

17. The following functions provide a standard interface to compiled code that has been linked into R: .C .Call

18. We will explore using .C and .Call with several code examples: Using .C I. Calling C with an integer vector II. Calling C with different vector types Using .Call III. Sending R integer vectors to C IV. Sending R character vectors to C V. Getting an integer vector from C VI. Getting a character vector from C VII. Getting a list from C

19. I. Calling C with an integer vector using .C

20. /* useC1.c */ void useC(int *i) { i[0] = 11; } The C function should be of type void. The compiled code should not return anything except through its arguments.

21. To compile the c code, type at the command prompt: R CMD SHLIB useC1.c The compiled code file name is useC1.so

22. In R: > dyn.load("useC1.so") > a <- 1:10 # integer vector > a [1] 1 2 3 4 5 6 7 8 9 10 > out <- .C("useC", b = as.integer(a)) > a [1] 1 2 3 4 5 6 7 8 9 10 > out$b [1] 11 2 3 4 5 6 7 8 9 10

23. You have to allocate memory to the vectors passed to .C in R by creating vectors of the right length. • The first argument to .C is a character string of the C function name. • The rest of the arguments are R objects to be passed to the C function.

24. All arguments should be coerced to the correct R storage mode to prevent mismatching of types that can lead to errors. • .C returns a list object. • The second .C argument is given the name b. This name is used for the respective component in the returned list object (but not passed to the compiled code).

25. II. Calling C with different vector types using .C

26. /* useC2.c */ void useC(int *i, double *d, char **c, int *l) { i[0] = 11; d[0] = 2.333; c[1] = "g"; l[0] = 0; }

27. To compile the c code, type at the command prompt: R CMD SHLIB useC2.c to get useC2.so To compile more than one c file: R CMD SHLIB file1.c file2.c file3.c to get file1.so

28. In R: > dyn.load("useC2.so") > i <- 1:10 # integer vector > d <- seq(length=3,from=1,to=2) # real number vector > c <- c("a", "b", "c") # string vector > l <- c("TRUE", "FALSE") # logical vector > i [1] 1 2 3 4 5 6 7 8 9 10 > d [1] 1.0 1.5 2.0 > c [1] "a" "b" "c" > l [1] "TRUE" "FALSE"

29. > out <- .C("useC", i1 = as.integer(a), d1 = as.numeric(d), c1 = as.character(c), l1 = as.logical(l)) > out $i1 [1] 11 2 3 4 5 6 7 8 9 10 $d1 [1] 2.333 1.500 2.000 $c1 [1] "a" "g" "c“ $l1 [1] FALSE FALSE

30. III. Sending R integer vectors to C using .Call

31. /* useCall1.c */ #include <R.h> #include <Rdefines.h> SEXP getInt(SEXP myint, SEXP myintVar) { int Imyint, n; // declare an integer variable int *Pmyint; // pointer to an integer vector PROTECT(myint = AS_INTEGER(myint));

32. Rdefines.h is somewhat more higher level then Rinternal.h, and is preferred if the code might be shared with S at any stage. • SEXP stands for Simple EXPression • myint is of type SEXP, which is a general type, hence coercion is needed to the right type. • R objects created in the C code have to be reported using the PROTECT macro on a pointer to the object. This tells R that the object is in use so it is not destroyed.

33. Imyint = INTEGER_POINTER(myint)[0]; Pmyint = INTEGER_POINTER(myint); n = INTEGER_VALUE(myintVar); printf(“ Printed from C: \n“); printf(“ Imyint: %d \n", Imyint); printf(“ n: %d \n", n); printf(“ Pmyint[0], Pmyint[1]: %d %d \n", Pmyint[0], Pmyint[1]); UNPROTECT(1); return(R_NilValue); }

34. The protection mechanism is stack-based, so UNPROTECT(n) unprotects the last n objects which were protected. The calls to PROTECT and UNPROTECT must balance when the user's code returns. • to work with real numbers, replace int with double and INTEGER with NUMERIC

35. In R: > dyn.load("useCall1.so") > myint<- c(1,2,3) > out<- .Call("getInt", myint, 5) Printed from C: Imyint: 1 n: 5 Pmyint[0], Pmyint[1]: 1 2 > out NULL

36. IV. Reading an R character vector from C using .Call

37. /* useCall2.c */ #include <R.h> #include <Rdefines.h> SEXP getChar(SEXP mychar) { char *Pmychar[5]; // array of 5 pointers // to character strings PROTECT(mychar = AS_CHARACTER(mychar));

38. // allocate memory: Pmychar[0] = R_alloc(strlen(CHAR(STRING_ELT(mychar, 0))), sizeof(char)); Pmychar[1] = R_alloc(strlen(CHAR(STRING_ELT(mychar, 1))), sizeof(char)); // ... and copy mychar to Pmychar: strcpy(Pmychar[0], CHAR(STRING_ELT(mychar, 0))); strcpy(Pmychar[1], CHAR(STRING_ELT(mychar, 1))); printf(“ Printed from C:”); printf(“ %s %s \n",Pmychar[0],Pmychar[1]); UNPROTECT(1); return(R_NilValue); }

39. In R: > dyn.load("useCall2.so") > mychar <- c("do","re","mi", "fa", "so") > out <- .Call("getChar", mychar) Printed from C: do re

40. V. Getting an integer vector from C using .Call

41. /* useCall3.c */ #include <R.h> #include <Rdefines.h> SEXP setInt() { SEXP myint; int *p_myint; int len = 5; // Allocating storage space: PROTECT(myint = NEW_INTEGER(len));

42. p_myint = INTEGER_POINTER(myint); p_myint[0] = 7; UNPROTECT(1); return myint; } // to work with real numbers, replace // int with double and INTEGER with NUMERIC

43. In R: > dyn.load("useCall3.so") > out<- .Call("setInt") > out [1] 7 0 0 0 0

44. VI. Getting a character vector from C using .Call

45. /* useCall4.c */ #include <R.h> #include <Rdefines.h> SEXP setChar() { SEXP mychar; PROTECT(mychar = allocVector(STRSXP, 5)); SET_STRING_ELT(mychar, 0, mkChar("A")); UNPROTECT(1); return mychar; }

46. In R: > dyn.load("useCall4.so") > out <- .Call("setChar") > out [1] "A" "" "" "" ""

47. VII. Getting a list from C using .Call

48. /* useCall5.c */ #include <R.h> #include <Rdefines.h> SEXP setList() { int *p_myint, i; double *p_double; SEXP mydouble, myint, list, list_names; char *names[2] = {"integer", "numeric"};

49. // creating an integer vector: PROTECT(myint = NEW_INTEGER(5)); p_myint = INTEGER_POINTER(myint); // ... and a vector of real numbers: PROTECT(mydouble = NEW_NUMERIC(5)); p_double = NUMERIC_POINTER(mydouble); for(i = 0; i < 5; i++) { p_double[i] = 1/(double)(i + 1); p_myint[i] = i + 1; }

50. // Creating a character string vector // of the "names" attribute of the // objects in out list: PROTECT(list_names = allocVector(STRSXP,2)); for(i = 0; i < 2; i++) SET_STRING_ELT(list_names,i,mkChar(names[i]));