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

[Douglas Bates]

An Rcpp object uses the storage allocated in R hence it has the same
column-major ordering of the elements in a matrix.  The same is true of
matrix objects in RcppArmadillo.

Basically anything that uses Lapack and BLAS or calling sequences
compatible with Lapack or BLAS uses column-major ordering.

With RcppEigen the default is column-major but you can specify row-major.
 However, when wrapping the Eigen object to return it to R there will be a
copying operation to convert row major to column major.

You are probably safest staying with column major ordering if possible.


On Wed, Nov 7, 2012 at 1:25 PM, Ramon Diaz-Uriarte <rdiaz02 at gmail.com>wrote:

>
> Dear All,
>
> I've recently started using Rcpp, so maybe this is a trivial question, but
> I have not been able to find an answer in previous posts or the docs.
>
>
> 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.
>
>
> --
> Ramon Diaz-Uriarte
> Department of Biochemistry, Lab B-25
> Facultad de Medicina
> Universidad Autónoma de Madrid
> Arzobispo Morcillo, 4
> 28029 Madrid
> Spain
>
> Phone: +34-91-497-2412
>
> Email: rdiaz02 at gmail.com
>        ramon.diaz at iib.uam.es
>
> http://ligarto.org/rdiaz
>
> _______________________________________________
> Rcpp-devel mailing list
> Rcpp-devel at lists.r-forge.r-project.org
> https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20121107/e01bf96f/attachment-0001.html>