[Rd] Nested SEXP functions

From: <statmobile_at_gmail.com>
Date: Thu, 15 Nov 2007 17:10:14 -0500

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


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_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Thu 15 Nov 2007 - 22:19:48 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 16 Nov 2007 - 02:30:18 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.