Frede Aakmann Tøgersen
2006-Mar-09 08:19 UTC
[R] When calling external C-function repeatedly I get differentresults; can't figure out why..
Not an expert in programming either, but to me it seems like you've forgotten to initialize the variable "tr". It just picks up garbage from allocated memory previously initialized by other processes. Med venlig hilsen Frede Aakmann T?gersen> -----Oprindelig meddelelse----- > Fra: r-help-bounces at stat.math.ethz.ch > [mailto:r-help-bounces at stat.math.ethz.ch] P? vegne af S?ren H?jsgaard > Sendt: 9. marts 2006 02:15 > Til: r-help at stat.math.ethz.ch > Emne: [R] When calling external C-function repeatedly I get > differentresults; can't figure out why.. > > 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 > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide! > http://www.R-project.org/posting-guide.html > >
Maybe Matching Threads
- When calling external C-function repeatedly I get different results; can't figure out why..
- question about dll crashing R
- Using matprod from array.c
- Speeding up matrix multiplies
- Internal state indicating if a data object has NAs/no NAs/not sure (Was: Re: Speeding up matrix multiplies)