[R] colSums in C

Peter Dalgaard BSA p.dalgaard at biostat.ku.dk
Sat Dec 15 01:08:59 CET 2001

David Brahm  <brahm at alum.mit.edu> writes:

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

Your code is not very efficiently written, but the compiler ought to
be able to catch the common subexpressions etc. On the other hand,
especially on Sparc there are some pretty tricky pipelining and
parallelism issues that can be really hard to get right and can mean
quite a lot. And there are some optimizations that a compiler cannot
really do without making assumptions that might be wrong.
(Prototypically: what if z and x have overlapping storage? *You* know
they don't, but the compiler really can't be sure.)
> 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];
> }

Hmm. What happens if you rewrite the "else" case as something like

register long i;
register double s, *p; 

p = x;
for(j=0; j<dx[1]; j++) {
        i = dx[0];
        s = 0.;
        while (i--)
                s += *p++;
        z[j] = s;

   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907
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

More information about the R-help mailing list