Array ordering in Julia

Question

Is there a way to work with C-ordered or non-contiguous arrays natively in Julia? For example, when using NumPy, C-ordered arrays are the default, but I can initialize a Fortran ordered array and do computations with that as well. One easy way to do this was to take the Transpose of a matrix. I can also work with non-contiguous arrays that are made via slicing. I have looked through the documentation, etc. and can't find a way to make, declare, or work with a C-ordered array in Julia. The transpose appears to return a copy.

Does Julia allow a user to work with C-ordered and non-contiguous arrays? Is there currently any way to get a transpose or a slice without taking a copy?

Edit: I have found how to do slicing. Currently it is available as a different type called a SubArray. As an example, I could do the following to get the first row of a 100x100 array A

sub(A, 1, 1:100)

It looks like there are plans to improve this, as can be seen in https://github.com/JuliaLang/julia/issues/5513

This still leaves open the question of C-ordered arrays. Is there an interface for C-ordered arrays? Is there a way to do a transpose via a view instead of a copy?

Answer 1

Naturally, there's nothing that prevents you from working with row-major arrays as a chunk of memory, and certain packages (like Images.jl) support arbitrary ordering of arbitrary-dimensional arrays.

Presumably the main issue you're wondering about is linear algebra. Currently I don't know of anything out-of-the-box, but note that matrix multiplication in Julia is implemented through a series of functions with names like A_mul_B, At_mul_B, Ac_mul_Bc, etc, where t means transpose and c means conjugate. The parser replaces expressions like A'*b with Ac_mul_B(A, b) without actually taking the transpose.

Consequently, you could implement a RowMajorMatrix <: AbstractArray type yourself, and set up special multiplication rules:

A_mul_B(A::RowMajorMatrix, B::RowMajorMatrix) = At_mul_Bt(A, B)
A_mul_B(A::RowMajorMatrix, B::AbstractArray) = At_mul_B(A, B)
A_mul_B(A::AbstractArray, B::RowMajorMatrix) = A_mul_Bt(A, B)

etc. In addition to these two-argument versions, there are 3-argument versions (like A_mul_B!) that store the result in a pre-allocated output; you'd need to implement those, too. Finally, you'd also have to set up appropriate show methods (to display them appropriately), size methods, etc.

Finally, Julia's transpose function has been implemented in a cache-friendly manner, so it's quite a lot faster than the naive

for j = 1:n, i = 1:m
    At[j,i] = A[i,j]
end

Consequently there are occasions where it's not worth worrying about creating custom implementations of algorithms, and you can just call transpose.

If you implement something like this, I'd encourage you to contribute it as a package, as it's likely that others may be interested.