[Rd] storage.mode, C data types and speed

[Heather Turner]

Hi,

I am trying to speed up part of an algorithm in which certain columns of a large matrix (X) are replaced by the element-wise product of a matrix (M) and a vector (v). In R, the code might be

X[, ind] <- M * v

I have written a small C routine to do this for me, but the timing depends on how I define the storage.mode of the objects in R and the data types in C, in a way which I don't understand.

To illustrate, suppose I have the following for X, M and v:
nr <- 10000
X <- matrix(1, nr, 500)
M <- matrix(1:(nr*450), nr, 450)
v <- 1:nr
storage.mode(X) <- storage.mode(M) <- storage.mode(v) <- "double"

Then I have the following integers required by my C code - these are the objects which I'm not sure how to handle
a <- 50*nr            #index of X indicating where to start replacing values
nm <- length(M)
nv <- length(v)

I would have thought I wanted

storage.mode(a) <- storage.mode(nm) <- storage.mode(nv) <- "integer"

to go with the following C code

# include <Rinternals.h>

SEXP prod_integer(SEXP X, SEXP M, SEXP v, SEXP a, SEXP nm, SEXP nv) {
  int i = INTEGER(a)[0], i1 = 0, i2 = 0;

  for ( ; i1 < INTEGER(nm)[0]; i2 = (++i2 == INTEGER(nv)[0]) ? 0 : i2) { 
    REAL(X)[i++] = REAL(M)[i1++] * REAL(v)[i2];
  }
  return(X);
}

Running this is R gives the following timings on my PC

> dyn.load("D:/C_routines/prod_integer")
> for(i in 1:3) {print(system.time(.Call("prod_integer", X, M, v, a, nm, nv)))}
[1] 0.17 0.00 0.18   NA   NA
[1] 0.18 0.00 0.17   NA   NA
[1] 0.15 0.00 0.17   NA   NA

But strangely (to me) if I change the storage mode of my integers to "double", I get

> storage.mode(a) <- storage.mode(nm) <- storage.mode(nv) <- "double"
> for(i in 1:3) {print(system.time(.Call("prod_integer", X, M, v, a, nm, nv)))}
[1]  0  0  0 NA NA
[1]  0  0  0 NA NA
[1]  0  0  0 NA NA

If on the other hand I use REAL instead of INTEGER in my C code:

# include <Rinternals.h>

SEXP prod_real(SEXP X, SEXP M, SEXP v, SEXP a, SEXP nm, SEXP nv) {
  int i = REAL(a)[0], i1 = 0, i2 = 0;

  for ( ; i1 < REAL(nm)[0]; i2 = (++i2 == REAL(nv)[0]) ? 0 : i2) { 
    REAL(X)[i++] = REAL(M)[i1++] * REAL(v)[i2];
  }
  return(X);
}

The reverse is true:
> storage.mode(a) <- storage.mode(nm) <- storage.mode(nv) <- "integer"
> for(i in 1:3) {print(system.time(.Call("prod_real", X, M, v, a, nm, nv)))}
[1]  0  0  0 NA NA
[1]  0  0  0 NA NA
[1]  0  0  0 NA NA
> storage.mode(a) <- storage.mode(nm) <- storage.mode(nv) <- "double"
> for(i in 1:3) {print(system.time(.Call("prod_real", X, M, v, a, nm, nv)))}
[1] 0.22 0.00 0.22   NA   NA
[1] 0.21 0.00 0.20   NA   NA
[1] 0.21 0.00 0.22   NA   NA
> identical(.Call("prod_integer", X, M, v, a, nm, nv), .Call("prod_real", X, M, v, a, nm, nv))
[1] TRUE

So I seem to get the fastest results if I store a, nm and nv as doubles in R and treat them as integers in C or if I store them as integers in R and treat them as doubles in C, rather than matching the storage in R to the data type in C.

I must be misunderstanding something here. Can someone explain what's going on - please note I have only just begun to learn C, apologies if this is a basic C issue.

Thanks,

Heather

Dr H Turner
Research Assistant
Dept. of Statistics
The University of Warwick
Coventry
CV4 7AL

Tel: 024 76575870
Url: www.warwick.ac.uk/go/heatherturner