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
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
