EIGEN_DEFAULT_TO_ROW_MAJOR and other flags

[HaoZeke]

I'm interfacing with a legacy C++ code which sets EIGEN_DEFAULT_TO_ROW_MAJOR (well known to break compatibility) but, basic searching and playing with defining the macro doesn't seem to help at all in R, matrices are still in column order and have to be manually reshaped (say matrix(byrow=T, ncol=3)).

I could wrap this thing in an R class and deal with this there every-time but was hoping I might have missed a way of setting the flag.

Things I tried:

#define EIGEN_DEFAULT_TO_ROW_MAJOR
#include <RcppEigen.h>
// [[Rcpp::depends(RcppEigen)]]

And moving it after the #include or the comment also does not help.

[eddelbuettel]

Can you disentangle this and maybe mock an example? Could it be that the EIGEN_DEFAULT_TO_ROW_MAJOR only affects representation in your C++ side whereas by the time we are back in R we always are an R matrix and hence bound by its representation? (That said I have seen 'have to transpose' situation eg when coming from NumPy over in another small repo of mine.

[HaoZeke]

Can you disentangle this and maybe mock an example? Could it be that the EIGEN_DEFAULT_TO_ROW_MAJOR only affects representation in your C++ side whereas by the time we are back in R we always are an R matrix and hence bound by its representation? (That said I have seen 'have to transpose' situation eg when coming from NumPy over in another small repo of mine.

Sorry, yes I should have thrown up a toy repo, will do so in a bit. The gist of the problem is how the matrix entries are filled with and without the flag. Consider the following:

...
            sscanf(line,"%lf %lf %lf\n", &x, &y, &z);
            positions(i,0) = x;
            positions(i,1) = y;
            positions(i,2) = z;
...

Which is meant to read simple things like coordinates:

14.776611 10.546847 10.871128
12.064038 11.226825 10.801342
15.192614 11.730865 10.909837
11.029134 11.970189 11.190553
13.079778 10.730201 11.722665
14.292513 12.611768 11.545665
10.264720 12.314561 10.501020
10.874573 12.279582 12.220608
13.061882 11.070527 12.764390 
14.178451 12.579462 12.629010
13.263563 9.648278 11.669731
14.106338 13.600457 11.139521
12.189980 10.941897 9.749161
16.147365 11.979995 10.436476

When the flag is not set or accounted for in R, Eigen fills entries in column order, so the matrix formed in R is:

           [,1]      [,2]      [,3]
 [1,] **14.776611** 11.722665 12.579462
 [2,] **10.546847** 14.292513 12.629010
 [3,] **10.871128** 12.611768 13.263563
 [4,] 12.064038 11.545665  9.648278
 [5,] 11.226825 10.264720 11.669731
 [6,] 10.801342 12.314561 14.106338
 [7,] 15.192614 10.501020 13.600457
 [8,] 11.730865 10.874573 11.139521
 [9,] 10.909837 12.279582 12.189980
[10,] 11.029134 12.220608 10.941897
[11,] 11.970189 13.061882  9.749161
[12,] 11.190553 11.070527 16.147365
[13,] 13.079778 12.764390 11.979995
[14,] 10.730201 14.178451 10.436476

When used with the flag internally on the C++ side, or when reshaped explicitly as mentioned with matrix(ncol=3, byRow=T) we get the correct data in R:

           [,1]      [,2]      [,3]
 [1,] 14.776611 10.546847 10.871128
 [2,] 12.064038 11.226825 10.801342
 [3,] 15.192614 11.730865 10.909837
 [4,] 11.029134 11.970189 11.190553
 [5,] 13.079778 10.730201 11.722665
 [6,] 14.292513 12.611768 11.545665
 [7,] 10.264720 12.314561 10.501020
 [8,] 10.874573 12.279582 12.220608
 [9,] 13.061882 11.070527 12.764390
[10,] 14.178451 12.579462 12.629010
[11,] 13.263563  9.648278 11.669731
[12,] 14.106338 13.600457 11.139521
[13,] 12.189980 10.941897  9.749161
[14,] 16.147365 11.979995 10.436476

P.S: Feel free to ignore this for now though, I'll link to a working runnable repo in a while.

P.P.S: I took a peak at the CNPy bindings (very neat) and yes, the effect of the Eigen flag is basically the same as the Fortran and C ordering in NDArrays (but set once per compile, instead of dynamically as is done in an NDArray)

[eddelbuettel]

Got it, thanks, this is now very clear.

But let me ask one more question though. Can you detail how / if our RcppEigen interact with your legacy software? I am thinking that if that part is built as it is (and independently from us over at R and RcppEigen) then we likely cannot change its behaviour -- but may have to rely to an appropriate tranpose() call. Does that make sense?