[Rcpp-devel] Problems returning fload

[Jeff Stevens]

Hi,

I'm trying to use Rcpp to speed up some nasty for loops in my R code
(Ubuntu 12.04, R 2.15.2, Rcpp 0.10.1, gcc 4.6.3).  However, I'm a newbie to
C and to Rcpp, and I'm having trouble returning the correct output.
Basically, I input an R matrix to an inline (using Rcpp plugin) C function
that runs through each element of the matrix and performs some calculations
to generate an array (float).  When I run the program in C, it calculates
the correct values.  However, when I run it in R, it give erratic results.
Here is example code that results in the same problem:
test_src <- '
Rcpp::NumericMatrix my_matrix(a);   // input matrix
int num_rows = my_matrix.nrow(), num_cols = my_matrix.ncol(); // extract
matrix dimensions
float my_sum[50];
float my_sum2[50];
Rcpp::NumericVector my_expr(num_cols);  // output array

    for(int i = 0; i < num_rows; i++) {
      for(int j = 0; j < num_cols; j++) {
      my_sum[j] = my_sum[j] + my_matrix(i,j);   // sum columns
      my_sum2[j] = (my_sum[j] + j);             // sum columns plus the
column index
      }
    }
    for(int j = 0; j < num_cols; j++) {
      my_expr[j] = my_sum[j] / my_sum2[j];      // divide sum by sum2
    }
return my_expr;
'
test_func <- cxxfunction(signature(a = "numeric"), test_src, plugin =
"Rcpp")

my_matrix <- matrix(c(0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 0), nrow =
5, ncol = 8)

test_func(my_matrix)

When I run this, I get:
1.0000000  0.6666667       NaN       NaN 1.0000000 0.2857143 0.3333333
0.3000000

If I run the following in C:
int num_rows = 5, num_cols = 8;
int my_matrix[5][8];
float my_sum[50];
float my_sum2[50];
float my_expr[50];
int my_matrix[5][8] = {
    {0, 1, 1, 0, 1, 0, 0, 1},
    {1, 0, 0, 0, 0, 0, 0, 1},
    {0, 0, 1, 0, 1, 1, 1, 1},
    {0, 1, 1, 0, 0, 0, 1, 1},
    {0, 0, 1, 0, 1, 1, 1, 0}
};

int main () {
    for(int i = 0; i < num_rows; i++) {
      for(int j = 0; j < num_cols; j++) {
      my_sum[j] = my_sum[j] + my_matrix[i][j];  // sum columns
      my_sum2[j] = (my_sum[j] + j);             // sum columns plus the
column index
      }
    }
    for(int j = 0; j < num_cols; j++) {
      my_expr[j] = my_sum[j] / my_sum2[j];      // divide sum by sum2
      printf("%f\t", my_expr[j]);
    }
}

I get
1.000000  0.666667   0.666667   0.000000   0.428571   0.285714   0.333333
0.363636

I'm sure that there is something simple that I am doing wrong.  Any
suggestions?

Thanks,
Jeff
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121219/b406bd84/attachment.html>