[R] question about dll crashing R

[Benn Fine]

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