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

[Ramon Diaz-Uriarte]

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