Support for row- and column-major order

[grothesque]

I reference a discussion from rust-ndarray/ndarray#1272 (comment):

@grothesque wrote:

Is your choice motivated by BLAS/LAPACK being (marginally) more efficient for column-major data?

Do I understand correctly that mdarray is column major in the sense that the restricted layouts are column major? But the fully strided layout can accept any (fixed rank) strided array, right? Right now in Rust we cannot have a fully generic ndspan like in C++, but it should be possible to have a set of useful layouts for both column-major and row-major within a single library, or do you see a problem with this?

@fre-hu replied:

The choice is only to have a convention, and then column major is common for linear algebra. It is used both for memory layout and to give the order of dimensions in iteration.

Using strided layout with row major data will work, but operations that depend on iteration order will have worse access pattern. It works fine for interfacing though, and internally one could make a copy or reverse indices.

To have full support for both row and column major would require one more generic parameter for the order. I had it in an earlier version, both removed it as it made both the library and interface more complex. C++ mdspan gets around this since it is quite thin.

From my point of view, row-major is arguably more relevant for a Rust array library than column major:

• Rust is much more a spiritual heir to C/C++ than Fortran.
• expr![[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]] having shape (3, 2) and not (2, 3) by default seems surprising.
• NumPy is and likely will remain the array library to which most people are first exposed. Despite its matlab lineage, but in consistency with the Python/C-world to which it belongs, NumPy uses row major by default. Moreover, seamless interoperability with NumPy (at least potentially) seems like an important feature of a Rust array library. (Cf. the success of Polars.)

So, if there can be only one, I'd vote for row major 😇...

However, as far as memory layout goes, it should be possible to have both without an additional generic parameter, right? Just like C++ mdspan has layout_right, layout_left, and layout_stride.

The problem seems to be more about ensuring efficient order of dimensions when iterating. One possibility would be to have both ("iterate_left_to_right", and "iterate_right_to_left"), and then only one (or none) would be efficient for a given array.

To treat the general case efficiently, there could be a function to (statically or dynamically) reorder dimensions into either layout (if possible).

All of this would not require an additional generic parameter (I believe).

[fre-hu]

I have switched to row-major order as it makes more sense, see the latest commit. I have also reduced the layout types to dense and strided, since these are the most important. Better to keep it simple and add back if it is really needed.

Regarding having both row and column-major order, yes the order could be merged into the layout type. The complexity is still there though, and there must be rules about iteration order, how to derive types for subarrays and how to do broadcasting etc. I will think some more about it.

[grothesque]

Sounds great. Looking forward to your ideas about merging order into layout.

I fully agree that an array crate for Rust should strive to be minimalist. I think that the right approach is to think of it as providing common concepts and glue, while actual numerical algorithms are to be provided by other crates.

This is the approach that has been successful with Fortran, and that C++ is finally pursuing with their mdspan. I consider the monolithic approach of Numpy/SciPy, despite their success, as a technical liability rooted in the shortcomings of Python's packaging.

[grothesque]

I have a question regarding the terms Dense and Strided (and also Flat and General that were removed). Are these terms standard or used in this way in other libraries?

I’m asking because Dense could be confusing as in numerics it often means the opposite of sparse. In this sense any mdarray array is dense.

[fre-hu]

Yes, Dense and Strided are from DenseArray and StridedArray in Julia. I have seen dense in more places to mean that elements are stored contiguously but I'm not sure if there is a clear definition. General comes from BLAS general matrix type.

[grothesque]

I believe that Julia’s use of the term “dense” does not reflect the usage in the wider community. Wouldn’t “contiguous” be the more appropriate term? (A web search for “dense matrix” vs “contiguous matrix” shows many examples for either term.) Alternatives could be `COrder` or `CContiguous`. This would keep the possibility for supporting Fortran-Order without using `Strided`. What triggered my above question is a student (who’s now using mdarray – thanks!) coming and asking whether mdarray has also sparse arrays, because he saw “dense”. So the potential for confusion is indeed there. My experience with teaching mdarray to this new user is that the recent renames help a lot with understanding things!

[fre-hu]

I see what you mean with dense vs sparse, and changing to contiguous could be an alternative. I will think some more on this, and perhaps together with other changes later on.

[grothesque]

Thanks. Meanwhile, I noticed that the term “contiguous” is already present in the API through the method is_contiguous.