Row Major Arrays

[tknopp]

I have already asked this on julia-users (https://groups.google.com/forum/?fromgroups=#!searchin/julia-users/Row$20Major/julia-users/FUU7C6Pvx3Q/Q6iUkUQM8xAJ) but maybe Github is the better place for this.

There are situations where row-major arrays are more suitable. In particular I am thinking about Kaczmarz algorithm which operates on matrix rows and this is simply more efficient on row-major arrays.

Solutions:

• Keep the reversed indices in mind when writing the algorithm. IMHO a hacky solution.
• Write some array view type which reverses indexing (i.e. the lazy transpose type)
• Let the standard array have different storage orders. This is how Numpy does it.

I know that this is core stuff of Julia and understand any objections in changing things here. Still, Julia seems to be not that flexible in this regards as Numpy and it would be great if we could get a clean solution for this.

[timholy]

Several existing algorithms follow LAPACK and allow you to pass 'N', 'T', or 'C' to indicate whether you want an operation to work with normal, transposed, or conjugate-transposed storage order. The algorithms take different code-paths for all 3 cases, and for efficiency I suspect one needs them all (cache-misses are so important that I think one needs to think about them explicitly).

For sparse matrices, this is an interesting approach, and perhaps the only viable one for making both A*v and A'*w fast.

[timholy]

CC @madeleineudell.

[tknopp]

I know that a flexible data layout would require a lot of work to get all operations efficient. But maybe one can use Julia fallbacks in the first iteration for row-major data.

[jiahao]

The Transpose{Array} type proposed in #42 could be a way of specifying a row-major matrix, although I'm not sure what people think of the corresponding CTranspose{Array} which would be needed to make use of Lapack 'C' storage.

[jiahao]

For sparse matrices, however, I see no conceptual reason to not support a SparseMatrixCSR type. Would be interesting to sketch out the possibilities for code reuse between CSR and CSC types.

[tknopp]

Yes a transposed array type would solve this although it is still somewhat inconvenient as one has to fill the array in a transposed way (My application has row major order matrices stored on disk). But having CSR in addition to CSC is the same as having column major + row major arrays.

[lindahua]

Having both CSR and CSC for sparse matrices are desirable.

However, for dense matrices, having two kinds of layouts would complicate things a lot:
(1) Element-wise operations that involve arrays of different layouts are very tricky to write.
(2) Lots of codes have been developed based on column-major order. Many of such codes won't work nicely or efficiently with row-major arrays.
(3) Shall we require/encourage people to consider both column-major or row-major cases when writing new functions that accept array arguments?

If we want to add the support of row-major arrays, it is much more preferable to add a RowMajorArray type than having Array to support both column-major and row-major orders.

[tknopp]

The RowMajorArray seems to be easiest approach to get something working and all your concerns are correct. I still wanted to put the flexible array layout solution onto the table as it is kind of more general.

I am wondering if I am the first to run into issues with row major data and Julia arrays. Numpy defaults to row major ordering.

[andreasnoack]

Isn't it mostly a matter about habit? If you are more used to row major you'd prefer to write your loop efficient for row major. Fortran is columns major and considered okay for numerical programming.

[lindahua]

I think Fortran and C respectively using column-major and row-major orders (by default) was unfortunate. Since then, languages are divided into two camps. According to Wikipedia:

Row-major order is used in C/C++, Mathematica, PL/I, Pascal, Python, Speakeasy, SAS and others. Column-major order is used in Fortran, OpenGL and OpenGL ES, MATLAB, GNU Octave, R, Julia, Rasdaman, and Scilab.

[tknopp]

@andreasnoackjensen: Technically this is not a problem. But my implementation of the Kaczmarz algorithm now operates on columns instead of rows. This makes the code less readable. Its not the biggest issue ever but whenever I see that Julia has a corner where it is not completely awesome I feel reponsible to report an issue on this :-)

[madeleineudell]

When this discussion is expanded to sparse arrays as well as dense ones, it becomes clear that there are algorithms for which one storage format is much faster than another. It's important to be able to support the needs of advanced users who care about these details while making sure that defaults are strong enough that casual users don't get confused.

For example, scipy supports 7 different sparse matrix formats, with very clear documentation on the strengths and weaknesses of each. So as a programmer, I can decide to eg build a matrix using coo format and convert it to csr for efficient A_b, or csc for efficient b'_A. But all these matrix types are clearly named as such, and actually raise warnings if you, say, try to set indices in a csr matrix.

[timholy]

But having CSR in addition to CSC is the same as having column major + row major arrays.

Not really. For dense arrays, memory locations are arranged in a predictable manner---calculating the address is a matter of arithmetic. For sparse arrays, it's completely different: there isn't any way to predict where something is stored without looking up the indexing. You're likely to get 2 cache misses just to find out where something is stored, and a third to fetch it.

[lindahua]

I wholeheartedly agree to support multiple sparse matrix formats. I personally have long been waiting for CSR in my own applications.

Having said that, the need of multiple layouts for dense arrays is not as obvious.

[andreasnoack]

I guess having more options is not always a good thing even though a more efficient solution exists. Left turns are inefficient in right-hand traffic but it is minor inefficiency compared to allowing mix of left and right-hand traffic.

The good news are that you have an efficient implementation for the transposed problem. Now you only have to write the version for the non-transposed problem:-)

The efficiency/simplicity tradeoff seems different for sparse, as pointed out by others. The only problem is that not even CSC is covered fully yet so there is a lot of work to be done.