You assume that when you call any of those functions its arguments are
protected (as it is the case when the function is called from within R).
However, I do not see why they would be if you do not do it explicitly.
Therefore, it is very likely that a result of one function, used in turn
in a call to the other one, is garbage collected on the way and what you
get can be totally wrong. The more R API calls you have within each of
those functions, the more likely it is that the unprotected argument is
garbage collected as the events are triggered within some of those
calls.
One should be very careful about protecting things in .Call interface.
Now here is what I mean in a coding example:
## R CODE #####
z1 = .Call("interfaceCFun1", as.numeric(x)) # OK
z2 = .Call("interfaceCFun12", as.numeric(x)) # Can be !OK
SEXP f1(SEXP);
// C code
SEXP interfaceCFun(SEXP x) {
/ /x is protected, so it prints smth and returns new var to R
return f1(x);
}
SEXP interfaceCFun2(SEXP x) {
// for first call to f1 it prints smth of x and returns new var
// to the second call, however in the second call the argument is no
// more x, it is no more protected and can be released by GC, so it is
// no clear what it is going to print there!
return f1(f1(x));
// this would be correct:
SEXP y;
PROTECT(y=f1(x));
// now when we print in the following function we have memory
protected
y = f1(y)
UNPROTECT(1);
return y;
}
SEXP f1(SEXP x) {
Rprintf("\n"); // can possibly trigger garbage collector
Rprintf(REAL(x)[0]);
return allocVector(REALSXP,2);
}
Best,
Oleg
On Thu, 2007-11-15 at 17:10 -0500, statmobile at gmail.com
wrote:> 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 at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel