R help - Aug 2006 - question about dll crashing R

I have ported some R code to C to make it faster.

I can perform .Call("foobar",....) once and it
works fine. Absolutely correct answer.

 If I put a loop inside foobar and run the  main code
routine more than 100 times, it crashes R.

Or if I call .Call("foobar"....) seperately more than
two tims it crashes R.


For the most part I am doing matirx multiplies using

EXP matrixprod (SEXP x , SEXP y )
{ int nrx , ncx , nry , ncy , mode;
SEXP xdims , ydims , ans ;
char *transa = "N" , *transb = "N" ;
double one = 1.0 , zero = 0.0 ;
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 ) ) ;
F77_CALL(dgemm) ( transa , transb , &nrx , &ncy , &ncx
, &one ,
REAL( x ) , &nrx , REAL( y ) , &nry , &zero ,
REAL( ans ) , &nrx ) ;
UNPROTECT( 1 ) ;
return( ans ) ;
}

I am also generating random multiavriate normals using
the (not pretty) code

*
 1) Generate P independent standard normal deviates -
Ei ~ N(0,1)
     2) Using Cholesky decomposition find A s.t.
trans(A)*A = COVM
     3) trans(A)E + MEANV ~ N(MEANV,COVM)
*/

SEXP mvntest (SEXP mean, SEXP cov, SEXP temp)


{ int nrx , ncx , nry , ncy ,info,mode;
SEXP xdims , ydims , ans;

int i,j, one=1;
info = 1;

xdims = getAttrib (mean , R_DimSymbol ) ;
ydims = getAttrib (cov , R_DimSymbol ) ;
mode = REALSXP;
nrx = INTEGER( xdims ) [ 0 ] ;
ncx = INTEGER( xdims ) [ 1 ] ;
nry = INTEGER( ydims ) [ 0 ] ;
ncy = INTEGER( ydims ) [ 1 ] ;



/* create the upper trianglular matrix A */

/* such that t(A) %*% A = Sigma */

GetRNGstate();	

F77_CALL(dpofa) ( REAL( cov ), &nry , &ncy , &info);
Rprintf("Info = %d\n",info);


for(i=0;i<nry;i++)
  for(j=0;j<i;j++)
    REAL(cov)[i+j*ncy] = 0.0;


PROTECT( ans = allocMatrix (mode, nrx , one ) ) ;
for(i=0;i<nry;i++)
  REAL(temp)[i] = rnorm(0,1);
ans = tmatrixprod(cov,temp);
for(i=0;i<nry;i++)
  REAL(ans)[i] = REAL(ans)[i]+REAL(mean)[i];
UNPROTECT( 1 ) ;
PutRNGstate();
return( ans ) ;


}


I have a feeling I am messing up memory usage
somewhere but haven't a clue. Do I need to do garbage
collecting inside the C program ?The fact that the
code
runs a few times before R crashes is driving me nuts. 
I send most of what I need into the C routine from R,
so I am not creating that many SEXP objects within the
program.

Any hints or ideas ?

Thanks!

Benn

On Thu, 3 Aug 2006, Benn Fine wrote:
> I have ported some R code to C to make it faster.
>
> I can perform .Call("foobar",....) once and it
> works fine. Absolutely correct answer.
>
> If I put a loop inside foobar and run the  main code
> routine more than 100 times, it crashes R.
>
> Or if I call .Call("foobar"....) seperately more than
> two tims it crashes R.
>
<snip>
> SEXP mvntest (SEXP mean, SEXP cov, SEXP temp)
>
>
> { int nrx , ncx , nry , ncy ,info,mode;
> SEXP xdims , ydims , ans;
>
> int i,j, one=1;
> info = 1;
>
> xdims = getAttrib (mean , R_DimSymbol ) ;
> ydims = getAttrib (cov , R_DimSymbol ) ;
> mode = REALSXP;
> nrx = INTEGER( xdims ) [ 0 ] ;
> ncx = INTEGER( xdims ) [ 1 ] ;
> nry = INTEGER( ydims ) [ 0 ] ;
> ncy = INTEGER( ydims ) [ 1 ] ;
>
>
>
> /* create the upper trianglular matrix A */
>
> /* such that t(A) %*% A = Sigma */
>
> GetRNGstate();
>
> F77_CALL(dpofa) ( REAL( cov ), &nry , &ncy , &info);
> Rprintf("Info = %d\n",info);
>
>
> for(i=0;i<nry;i++)
>  for(j=0;j<i;j++)
>    REAL(cov)[i+j*ncy] = 0.0;
>
>
> PROTECT( ans = allocMatrix (mode, nrx , one ) ) ;
> for(i=0;i<nry;i++)
>  REAL(temp)[i] = rnorm(0,1);
> ans = tmatrixprod(cov,temp);
^^^^^^^^^
Here you are returning a new SEXP from tmatrixprod but not protecting it. 
Remember, PROTECT() operates on the *value* of its argument, so it 
protects the thing that ans points to. When ans points to a new thing, it 
is still the old thing that is protected.

 	-thomas
> for(i=0;i<nry;i++)
>  REAL(ans)[i] = REAL(ans)[i]+REAL(mean)[i];
> UNPROTECT( 1 ) ;
> PutRNGstate();
> return( ans ) ;
>
>
> }
>
>
> I have a feeling I am messing up memory usage
> somewhere but haven't a clue. Do I need to do garbage
> collecting inside the C program ?The fact that the
> code
> runs a few times before R crashes is driving me nuts.
> I send most of what I need into the C routine from R,
> so I am not creating that many SEXP objects within the
> program.
>
> Any hints or ideas ?
>
> Thanks!
>
> Benn
>
> ______________________________________________
> 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
> and provide commented, minimal, self-contained, reproducible code.
>
Thomas Lumley			Assoc. Professor, Biostatistics
tlumley at u.washington.edu	University of Washington, Seattle