C code called by R continues to fail

Question

The following is part of the C code I wrote. The foo function should be called in R. The code continues to call R, and I narrowed down the problem to this outer() function, which is used to calculate the external sum or difference. Pay attention to the part that is commented out: if I do not comment on this, the function will cause R to crash if each of the arrays contains, say, more than 1000 data points . If I comment on this, I can calculate the external sum / difference for significantly longer arrays without problems (for example, more than 100,000 data points per array). I wonder what the problem is ... Thanks!

 #include <Rh> #include <Rmath.h> #include <stdio.h> #include <math.h> #include <stdlib.h> void outer(double *x1, double *x2, int *n, int operation, double *output){ int i, j; if(operation==1){ for(i=0; i<*n; i++){ for(j=0; j<*n; j++){ output[(*n)*i+j]=x1[j]+x2[i]; } } } else if(operation==2){ for(i=0; i<*n; i++){ for(j=0; j<*n; j++){ output[(*n)*i+j]=x1[j]-x2[i]; //Rprintf("%d ", (*n)*i+j); //<-----------HERE } } } } void foo(double *x, double *y, int *npred, int *nsamp){ int oper=2; double xouter[*nsamp], youter[*nsamp]; double outer_temp_x[(*nsamp)*(*nsamp)], outer_temp_y[(*nsamp)*(*nsamp)]; outer(x, x, nsamp, oper, &outer_temp_x[0]); outer(y, y, nsamp, oper, &outer_temp_y[0]); } 

// After compiling the code, I use the code below in R to call the function:

 dyn.load("foo.so") x=as.matrix(rnorm(10000)) y=rlnorm(10000) invisible(.C("foo", x=as.double(as.vector(x)), y=as.double(y), npred=as.integer(ncol(x)), nsamp=as.integer(length(y)) ) 

Answer 1

I think this is superior to the stack and creates problems.

Try the following:

 void foo(double *x, double *y, int *npred, int *nsamp){ int oper=2; double xouter[*nsamp], youter[*nsamp]; // The prior code allocated on the stack. Here, we make a pair of calls // to 'malloc' to allocate memory for the arrays. This gets memory from // the heap. The stack is fairly limited, but the heap is huge. // 'malloc' returns a pointer to the allocated memory. double* outer_temp_x=malloc(sizeof(double)*(*nsamp)*(*nsamp)); double* outer_temp_y=malloc(sizeof(double)*(*nsamp)*(*nsamp)); outer(x, x, nsamp, oper, &outer_temp_x[0]); outer(y, y, nsamp, oper, &outer_temp_y[0]); // The downside of allocating on the heap, is that you must release the // memory at some point. Otherwise you have what called a "memory leak." // 'free' is the function to free the memory, and it is called on the // pointer value returned by 'malloc'. free(outer_temp_x); free(outer_temp_y); }