[R] colSums in C

[David Brahm]

I wrote:
> I asked how to write speedy C code to implement colSums()...
> I've taken Doug [Bates]'s code, added names to the result, and included an
> na.rm flag.  Unfortunately, my na.rm option makes it really slow again
> ... no faster than pre-processing the matrix with "m[is.na(m)] <- 0".

Well, I'm an idiot.  I had forgotten to remove the R code that does that
pre-processing, so no wonder it wasn't any faster!  After fixing that, I got
results closer to those reported by Thomas Lumley, about 2.18s with ISNA, and
even faster (1.41s) with isnan (Peter Dalgaard's suggestion).  Is there any
reason to distrust isnan()?

Current colSums code is included below, and there's a rowSums too.  I'm still
prettying it up, but I'd like to open discussion on what to do with it when
it's ready for prime time:
  1) I can just use it locally and be perfectly happy,
  2) I can upload to CRAN (with Doug's permission, of course),
  3) I believe Doug Bates and Kurt Hornik once discussed modifying apply() to
     catch these special cases and use the optimized code.  Dangerous?

Thanks, everyone!

---- begin "colSums.R" ----
colSums <- function(x, na.rm=F) .Call("colSums", x, na.rm)
---- end "colSums.R" ----

---- begin "colSums.c" (90% due to Doug Bates) ----
#include <R.h>
#include <Rdefines.h>

SEXP colSums(SEXP m, SEXP narm) {
  register int i, j, NaRm;
  int *mdims, n, p;
  double *mm, sum;
  SEXP val, nms;

  if (!isMatrix(m)) error("m must be a matrix");
  mdims = INTEGER(GET_DIM(m));
  n = mdims[0]; p = mdims[1];
  PROTECT(val  = NEW_NUMERIC(p));
  PROTECT(nms  = GET_COLNAMES(GET_DIMNAMES(m)));
  PROTECT(  m  = AS_NUMERIC(m));
  NaRm = asLogical(narm);
  mm = REAL(m);
  for (j = 0; j < p; j++) {
    for (sum = 0., i = 0; i < n; i++) if (!NaRm || !isnan(mm[i])) sum += mm[i];
    REAL(val)[j] = sum;
    mm += n;
  }
  if (!isNull(nms)) namesgets(val, nms);
  UNPROTECT(3);
  return val;
}
---- end "colSums.c" ----

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