[Rcpp-devel] Matrix: column-major storage?

[Romain Francois]

Le 07/11/12 20:25, Ramon Diaz-Uriarte a écrit :
>
> Dear All,
>
> I've recently started using Rcpp

welcome.

> so maybe this is a trivial question, but
> I have not been able to find an answer in previous posts or the docs.

Rcpp plays with the rules set by R. And R says : column major.

> I thought that, for any of the Matrix types available, if we have
>
> x(i, j)
>
> then
>
> x(i, j + 1)
>
> would be contiguous in memory (whereas x(i + 1, j) could be extremely far
> away, depending on the size of the matrix).
>
> However, it seems that with Matrix objects it is the other way around, so
> that x(i, j) is next to x(i + 1, j), not to x(i, j + 1):
>
> #################
>
> showThePointers <- cxxfunction(signature(dm = "integer"),
>                       body = '
>
> int m = as<int>(dm);
> IntegerMatrix xm(m, m);
> int *p;
> int i;
> int j;
>
> std::cout << std::endl << "***** Row-major access *****" << std::endl;
>
> for(i = 0; i < m; i++) {
>     for(j = 0; j < m; j++) {
>       xm(i, j) = i * m + j;
>       p = &xm(i, j);
>       std::cout << "i = " << i << "; j = " <<
>       j << "; xm(i, j) = " << xm(i, j) << "; p " << p << std::endl;
>    }
> }
>
> std::cout << std::endl << " ****** Column-major access *****" << std::endl;
>
> for(j = 0; j < m; j++) {
>     for(i = 0; i < m; i++) {
>       xm(i, j) = i * m + j;
>       p = &xm(i, j);
>       std::cout << "i = " << i << "; j = " <<
>       j << "; xm(i, j) = " << xm(i, j) << "; p " << p << std::endl;
>    }
> }
>
> ',
>                       plugin = "Rcpp",
>                       verbose = TRUE)
>
> showThePointers(3)
>
> #################
>
>
>
>
> I noticed that because if we time functions that access and modify in
> row-major order, they are a lot slower than in column-major order (in my
> laptop, about four times slower in the examples below, about nine times if
> I then return the object to R).
>
>
> #################
>
> byRow <- cxxfunction(signature(dm = "integer"),
>                       body = '
>
> int m = as<int>(dm);
> IntegerMatrix xm(m, m);
> int i;
> int j;
> for(i = 0; i < m; i++) {
>     for(j = 0; j < m; j++) {
>       xm(i, j) = i * m + j;
>    }
> }
> ',
>                       plugin = "Rcpp",
>                       verbose = TRUE)
>
> byCol <- cxxfunction(signature(dm = "integer"),
>                       body = '
>
> int m = as<int>(dm);
> IntegerMatrix xm(m, m);
> int i;
> int j;
> for(j = 0; j < m; j++) {
>     for(i = 0; i < m; i++) {
>       xm(i, j) = i * m + j;
>    }
> }
> ',
>                       plugin = "Rcpp",
>                       verbose = TRUE)
>
>
> nn <- 10000
> benchmark(byRow(nn), byCol(nn),
>            columns=c("test", "replications", "elapsed", "relative",
>            "user.self", "sys.self"),
>            order="relative", replications=10)
>
>
> #################
>
>
> I understand it is not possible to store things in row-major order with
> Rcpp Matrix types?
>
>
> Best,
>
>
> R.
>
>


-- 
Romain Francois
Professional R Enthusiast
+33(0) 6 28 91 30 30

R Graph Gallery: http://gallery.r-enthusiasts.com
`- http://bit.ly/SweN1Z : SuperStorm Sandy

blog:            http://romainfrancois.blog.free.fr
|- http://bit.ly/RE6sYH : OOP with Rcpp modules
`- http://bit.ly/Thw7IK : Rcpp modules more flexible