[Rd] Nested SEXP functions

[statmobile at gmail.com]

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!