On 15 March 2017 at 14:17, Roland Denis wrote:
| Hello,
|
| while comparing the performance between pure R and different other
| syntaxes that mix R and Rcpp, I found out a strange behavior when using
| Rcpp sugar on NumericMatrix. As far as I understand Rcpp (I'm a newbie
| in R and Rcpp), syntactic sugar is designed for NumericVector, but since
| NumericMatrix seems to be seen like an augmented NumericVector, I've
| supposed this works on NumericMatrix too (and, in fact, it compiles).
|
| Here is a minimal example :
|
| library(Rcpp)
| cppFunction("
| NumericMatrix scaleSugar( NumericMatrix A, double k )
| {
| return k*A;
| }
| ")
|
| A <- diag(3)
| B = scaleSugar(A, 2)
|
| Displaying A and B returns :
| > A
| [,1] [,2] [,3]
| [1,] 2 0 0
| [2,] 0 2 0
| [3,] 0 0 2
| > B
| [,1] [,2] [,3]
| [1,] 2 0 0
| [2,] 0 2 0
| [3,] 0 0 2
|
| As you see, the matrix A seems to be multiplied in-place and returned.
|
| Creating the returned matrix before the multiplication doesn't change
| the result :
| library(Rcpp)
| cppFunction("
| NumericMatrix scaleSugar( NumericMatrix A, double k )
| {
| NumericMatrix R( A.nrow(), A.ncol());
| R = k*A;
| return R;
| }
| ")
|
| A <- diag(3)
| B = scaleSugar(A, 2)
|
| > A
| [,1] [,2] [,3]
| [1,] 2 0 0
| [2,] 0 2 0
| [3,] 0 0 2
|
| However, copying A into R before the scaling works :
| library(Rcpp)
| cppFunction("
| NumericMatrix scaleSugar( NumericMatrix A, double k )
| {
| NumericMatrix R( A.nrow(), A.ncol(), A.begin());
| return k*R;
| }
| ")
|
| A <- diag(3)
| B = scaleSugar(A, 2)
|
| > A
| [,1] [,2] [,3]
| [1,] 1 0 0
| [2,] 0 1 0
| [3,] 0 0 1
|
|
| and, *more surprisingly*, reversing the multiplication operand's order
| works too :
| library(Rcpp)
| cppFunction("
| NumericMatrix scaleSugar( NumericMatrix A, double k )
| {
| return A*k;
| }
| ")
|
| A <- diag(3)
| B = scaleSugar(A, 2)
|
| > A
| [,1] [,2] [,3]
| [1,] 1 0 0
| [2,] 0 1 0
| [3,] 0 0 1
|
| Are Rcpp sugars supposed to work with NumericMatrix ?
That could be a bug. But please note that we said for many years that the
matrix _calculations_ are underdeveloped / incomplete (see eg the Rcpp FAQ
suggesting to do all this with RcppArmadillo and Armadillo matrices). So if
I were you I'd take a good look at RcppArmadillo.
Dirk
--
http://dirk.eddelbuettel.com | @eddelbuettel | e...@debian.org
_______________________________________________
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
