Hey All, I was wondering if I could solicit a little advice. I have been experiencing some quirkiness in my C code through .Call. Unfortunately, my program is rather large, so I'm trying to create a brief example. Yet before I do, I thought maybe a conceptual question would be sufficient.
Basically, I'm writing up a multidimensional Gaussian likelihood (with spatial covariances and other enhancements). What I'm noticing is that when I have nested SEXP functions I get inconsistent results. Here is what I am essentially trying to accomplish when looking at the Gaussian kernel: l(beta) = (y-X*beta)^T V^{-1} (y-X*beta) Now in order to accomplish this, I wrote various linear algebra subroutines to handle the R objects, we'll call this: SEXP XtimesY(SEXP X,SEXP Y); // X*Y SEXP XminusY(SEXP X,SEXP Y); // X-Y SEXP tX(SEXP X); // X^T SEXP mycholinv(SEXP V); // Use cholesky decomposition for inverse Now, what I'm noticing is that if I call each routine individually such as: pt1=XtimesY(X,beta); // X*beta pt2=Xminus(Y,pt1); // Y-X*beta pt3=tX(pt2); // (Y-X*beta)^T pt4=mycholinv(V); //V^{-1} pt5=XtimesY(pt2,pt4); // (Y-X*beta)^T V^{-1} result=XtimesY(pt5,pt2); //(y-X*beta)^T V^{-1} (y-X*beta) Then the result is correct. But if instead I want to save some lines of code, and I use: result=XtimesY(XtimesY(tX(XminusY(Y,XtimesY(X,beta))),mycholinv(V)),XminusY(Y,XtimesY(X,beta))) I don't always get the same result. Now my question is, should I expect weird and ugly things to happen when nesting SEXP functions such as this? Or is it even highly discouraged in general in C to do something like this? If this should work, then I'll need to go through each one of the functions and try to diagnose where the problem lies. Yet if it shouldn't work, I'll stick with the first approach I have going now. Thanks in advance for your input! ______________________________________________ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel