[R] When calling external C-function repeatedly I get different results; can't figure out why..

[Søren Højsgaard]

Dear all,
I need to calculate tr(xyxy) fast for matrices x and y. Inspired by the R-source code, I've created the following functions (I am *new* to writing external C-functions, so feel free to laugh at my code - or, perhaps, suggest changes):
 
#include <Rinternals.h>
#include <R_ext/Applic.h> /* for dgemm */

static void matprod(double *x, int nrx, int ncx,
      double *y, int nry, int ncy, double *z)
{
    char *transa = "N", *transb = "N";
    double one = 1.0, zero = 0.0;
    F77_CALL(dgemm)(transa, transb, &nrx, &ncy, &ncx, &one,
      x, &nrx, y, &nry, &zero, z, &nrx);
}

SEXP trProd2(SEXP x, SEXP y)
{
  int nrx, ncx, nry, ncy, mode, i;
  SEXP xdims, ydims, ans, ans2, tr;
  xdims = getAttrib(x, R_DimSymbol);
  ydims = getAttrib(y, R_DimSymbol);
  mode = REALSXP;
  nrx = INTEGER(xdims)[0];
  ncx = INTEGER(xdims)[1];
  nry = INTEGER(ydims)[0];
  ncy = INTEGER(ydims)[1];
  PROTECT(ans  = allocMatrix(mode, nrx, ncy));
  PROTECT(ans2 = allocMatrix(mode, nrx, ncy));
  PROTECT(tr   = allocVector(mode, 1));
  matprod(REAL(x), nrx, ncx, REAL(y), nry, ncy, REAL(ans)); 
  
  matprod(REAL(ans), nrx, ncy, REAL(ans), nrx, ncy, REAL(ans2)); 
  
  for (i=0; i< nrx; i++){
    REAL(tr)[0] = REAL(tr)[0] + REAL(ans2)[i*(nrx+1)];
  }
  UNPROTECT(3);
  return(tr);
}

In R I do:
 
trProd2 <- function(x, y) {
.Call("trProd2", x, y)
}
x <- y <- matrix(1:4,ncol=2)
storage.mode(x) <- storage.mode(y) <- "double"

for (i in 1:10)
  print(trProd2(x, y))
[1] 835
[1] 833
[1] 833
[1] 862
[1] 834
[1] 835
[1] 834
[1] 836
[1] 862
[1] 833

The correct answer is 833. Can anyone give me a hint to what is wrong? I am on a windows xp platform.
Thanks in advance
Søren