[R] colSums in C

[David Brahm]

Hi, all!

   My project today is to write a speedy colSums(), which is a function
available in S-Plus to add the columns of a matrix.  Here are 4 ways to do it,
with the time it took (elapsed, best of 3 trials) in both R and S-Plus:

m <- matrix(1, 400, 40000)
x1 <- apply(m, 2, sum)                  ## R=16.55   S=52.39
x2 <- as.vector(rep(1,nrow(m)) %*% m)   ## R= 2.39   S= 8.52
x3 <- g.apply.sum(m,2,F)                ## R= 5.72   S=11.25
x4 <- colSums(m,F)                      ##           S= 0.48

Method 1 is the naive "apply", and is significantly faster in R than S-Plus.
Method 2 is a matrix operation I learned from Bill Venables; also faster in R.
Method 3 is my C code, given below.
Method 4 is the S-Plus built-in function, the bogey to beat!

My question is, can anyone suggest how to speed up my C code to get results
comparable to colSums, or even to method 2?  Is it stupid code, a bad compiler,
non-optimal optimizer options, or ???

Here's my interface function and simple C code (compiled with gcc via R
SHLIB, no Makevars file):

g.apply.sum <- function(x, MARGIN=1, na.rm=T) {
  dx <- dim(x)
  if (na.rm) x[is.na(x)] <- 0
  .C("applysum", as.real(x), dx, as.integer(MARGIN), z=real(dx[MARGIN]),
     DUP=FALSE)$z
}

void applysum(double *x, long *dx, long *mar, double *z) {
  long i, j;
  if (*mar==1)
    for(i=0; i<dx[0]; i++) for(j=0; j<dx[1]; j++) z[i] += x[i + dx[0]*j];
  else
    for(j=0; j<dx[1]; j++) for(i=0; i<dx[0]; i++) z[j] += x[i + dx[0]*j];
}

Possibly useful audit trail:

agate|Rpkg> R SHLIB g.colSums/src/*.c
/res/local/bin/gcc -I/res/local/lib/R/include  -I/res/local/include   -fPIC  -g -O2 -c g.colSums/src/applyfilt.c -o g.colSums/src/applyfilt.o
/res/local/bin/gcc -I/res/local/lib/R/include  -I/res/local/include   -fPIC  -g -O2 -c g.colSums/src/applysum.c -o g.colSums/src/applysum.o
/res/local/bin/gcc -shared  -o g.colSums/src/applyfilt.so g.colSums/src/applyfilt.o g.colSums/src/applysum.o  -L/res/local/include  

Version:
 platform = sparc-sun-solaris2.6
 arch = sparc
 os = solaris2.6
 system = sparc, solaris2.6
 status = 
 major = 1
 minor = 3.1
 year = 2001
 month = 08
 day = 31
 language = R

-- 
                              -- David Brahm (brahm at alum.mit.edu)
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._