JuliaLang/julia — pull request #8432 by timholy (2014-09-21 01:02:56+00:00)
Title: Efficient cartesian iteration (new version of #6437)

This is the complement of "adding array views": while it seems not to be widely appreciated, the reality is that none of the View proposals (mine included) have good performance in the general case _if_ we rely on linear indexing (demos below, see for example also a [recent StackOverflow post](http://stackoverflow.com/questions/25872585/alternative-to-scale-for-subarray-objects)). This is basically a direct consequence of the fact that `div` is horribly slow. Consequently, we need to transition to writing algorithms with multidimensional/cartesian iterators. While the Cartesian macros provide a way to do this, I think it's fair to say that everyone would like a solution that doesn't rely on macros. While this PR offers only a small fraction of the power of Cartesian, at least simple cases can now get good performance without resorting to macros.

This is basically a redo of the long-languishing #6437, with two differences:
- It exploits the `@inline` macro I just merged
- It's based on @carlobaldassi's [`immutable` approach](https://gist.github.com/carlobaldassi/9413784) rather than the tuple approach I tried in #6437. With the tuples there was a 2.5x penalty, which I think is enough in practice to discourage use.

At least for low dimensional `Array`s, which have efficient linear indexing, there's scarcely any difference between linear and cartesian iteration. Nevertheless I preserved linear indexing for `Array`s using the ["traits trick"](https://github.com/JuliaLang/julia/issues/2345#issuecomment-54537633). If you have an `AbstractArray` that has fast linear indexing, you can cause `start` to use it simply by declaring `Base.linearindexing(::MyArrayType) = Base.LinearFast()`.

The only user-visible disadvantage I'm aware of in not using tuples is that the notation for multidimensional indexing needs to be `A[I]` rather than `A[I...]`. Consequently if/when we switch to tuples, people will have to update their code. In my opinion, that's a small price to pay for finally having efficient multidimensional iteration.

In the performance demos below, there's still room for some improvement. The most important one, I suspect, would be for types that have tuples in them to "inline" the tuple into the type (as pointed out by @simonster among others). Of course if we get that, this could be redone using tuples.

Performance demos:
```julia
function mdsum(A)
s = 0.0
for a in eachelement(A)
s += a
end
s
end

function mdsum2(A)
s = 0.0
@inbounds for I in eachindex(A)
s += A[I]
end
s
end

function sumlinear(A)
s = 0.0
@inbounds for i = 1:length(A)
s += A[i]
end
s
end

function sumcart2(A::AbstractMatrix)
s = 0.0
@inbounds for j = 1:size(A,2), i = 1:size(A,1)
s += A[i,j]
end
s
end

println("Arrays")
A = rand(10^4, 10^4)
mdsum(A); mdsum2(A); sumlinear(A); sumcart2(A)
@time mdsum(A)
@time mdsum2(A)
@time sumlinear(A)
@time sumcart2(A)

println("SubArrays")
B = sub(A, 2:9999, 2:9999)
mdsum(B); mdsum2(B); sumlinear(B); sumcart2(B)
@time mdsum(B)
@time mdsum2(B)
@time sumlinear(B)
@time sumcart2(B)

using ArrayViews
println("ArrayViews")
B = view(A, 2:9999, 2:9999)
mdsum(B); mdsum2(B); sumlinear(B); sumcart2(B)
@time mdsum(B)
@time mdsum2(B)
@time sumlinear(B)
@time sumcart2(B)
```
Results:
```julia
julia> include("/tmp/testiteration.jl")
Arrays
elapsed time: 0.15763366 seconds (13896 bytes allocated)
elapsed time: 0.209725186 seconds (96 bytes allocated)
elapsed time: 0.133438424 seconds (96 bytes allocated)
elapsed time: 0.133621887 seconds (96 bytes allocated)
SubArrays
elapsed time: 0.491196255 seconds (96 bytes allocated)
elapsed time: 0.406468852 seconds (96 bytes allocated)
elapsed time: 1.871101162 seconds (96 bytes allocated)
elapsed time: 0.283826772 seconds (96 bytes allocated)
ArrayViews
elapsed time: 1.331280489 seconds (96 bytes allocated)
elapsed time: 0.27865117 seconds (96 bytes allocated)
elapsed time: 2.543739468 seconds (96 bytes allocated)
elapsed time: 0.253180521 seconds (96 bytes allocated)
4.99765960924932e7
```