[Rcpp-devel] Strange behaviro using Rcpp sugar with NumericMatrix

[Roland Denis]

Ok, thanks, I will take a look at the Armardillo and Eigen wrappers !

-- 
Roland


Le 15/03/2017 à 15:13, Dirk Eddelbuettel a écrit :
> 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
>
>
>