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

[Frede Aakmann Tøgersen]

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