for j = 1:3, i = 1:2
stuff with A[i,j]
end
Yes, this is an unfortunate consequence of mathematics being column-major – oh how I wish it weren't so.
On Thursday, October 29, 2015 at 1:24:23 PM UTC-4, Stefan Karpinski wrote:Yes, this is an unfortunate consequence of mathematics being column-major – oh how I wish it weren't so.I'm not sure what you mean by "mathematics being column major". Fortran is column-major, and a lot of the linear-algebra libraries (LAPACK etc.) are column-major in consequence, and Matlab is column-major.However, if Julia had chosen row-major, it would have been possible with a little care to call LAPACK and other column-major libraries without making a copy of any data, basically because for every linear-algebra operation there is an equivalent operation on the transpose of the matrix.
Yes, this is an unfortunate consequence of mathematics being column-major – oh how I wish it weren't so. The storage order is actually largely irrelevant, the whole issue stems from the fact that the element in the ith row and the jth column of a matrix is indexes as A[i,j]. If it were A[j,i] then these would agree (and many things would be simpler). I like your explanation of "an index closer to the expression to be evaluated runs faster" – that's a really good way to remember/explain it.
in my experience practical applications favor row ordering
The keyword there is that it was your experience. The thing about something as central to number crunching as organizing memory layouts in row- or column major orderings, is that there are more practical applications than any of us can imagine - in fact, most likely we couldn’t come up with all possible cases where one is significantly better than the other even if we stopped discussing Julia and just discussed that for a year.
For some applications column major is better, but for others row major rocks. We have to choose, and even if the choice was once quite arbitrary, the gain from switching is so extremely low it’s ridiculous to even consider it. Thus, I don’t think it’s fruitful to discuss here.
However, Julia is an incredibly expressive language, and you could quite easily write a macro that re-orders the indices for you:
julia> macro rowmajor(expr)
@assert expr.head == :ref
Expr(:ref, expr.args[1], expr.args[end:-1:2]...)
end
julia> macroexpand(:(@rowmajor A[i,j]))
:(A[j,i])
Of course, this only works on one indexing expression at a time, but it can be extended to work recursively through an expression tree and rewrite all array indexing expressions it finds, allowing usage as
@rowmajor begin
# formulate row-major algorithm here
end
No need to change the language - just bend it to do what you need :)
// T
We have to choose
for i2=1:n,
i1=1:m
A[ i1, i2 ] = 10 * i1 + i2
end
A = [ 10 * i1 + i2 for i2=1:n, i1=1:m ]