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
dyn.load("r_gsl_test.so")
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)
max(abs(c-d))

> 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.

### Like this:

Like Loading...

*Related*