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
_______________________________________________
Rcpp-devel mailing list
Rcpp-devel@lists.r-forge.r-project.org
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
