Linking R to GSL

I put the following in r_gsl_test.c

#include <R.h>
#include <stdio.h>
#include <gsl/gsl_blas.h>
void mmult(double *a, double *b, int *arow, int *acol, int *bcol, double *c) {
   gsl_matrix_view A = gsl_matrix_view_array(a, *arow, *acol);
   gsl_matrix_view B = gsl_matrix_view_array(b, *acol, *bcol);
   gsl_matrix_view C = gsl_matrix_view_array(c, *arow, *bcol);
   /* Compute C = A B */
   gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                       1.0, &A.matrix, &B.matrix,
                       0.0, &C.matrix); 

Then ran PKG_LIBS="-lgsl" R CMD SHLIB r_gsl_test.c at the command line (obviously, I need to have GSL installed).

Then tested using the following R code:

# PKG_LIBS="-lgsl" R CMD SHLIB r_gsl_test.c

mmult <- function(a, b) {
  z <- .C("mmult", as.double(t(a)), as.double(t(b)), nrow(a),
                    ncol(a), ncol(b), results=double(nrow(a)*ncol(b)))
  matrix(z$results, nrow=nrow(a), byrow=TRUE)
n <- 1000; k <- 100; m <- 200

a <- matrix(rnorm(k*m), nrow=k)
b <- matrix(rnorm(m*n), nrow=m)
system.time(c <- mmult(a,b))
system.time(d <- a %*% b)
> system.time(c <- mmult(a,b))
   user  system elapsed 
  0.016   0.002   0.037 
> system.time(d <- a %*% b)
   user  system elapsed 
  0.027   0.000   0.074 
> max(abs(c-d))
[1] 6.394885e-14

No huge speed advantage evident here, though it is proof of the concept of linking R to GSL.

This entry was posted in R. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s