JuliaLang/julia — pull request #3100 by toivoh (2013-05-14 18:26:29+00:00)
Title: RFC, WIP: Staged bsxfun and other broadcast operations

This patch adds an implementation of `bsxfun(op::Function, As::Array...)`
that generates inner loop functions depending on `op`, the number of arguments
`As`, and the dimension of the result.

The machinery is also used to provide
* broadcasting operators `.+ .- .* ./` for `Array` arguments
* `bsxfun!(op::Function, As::Array...)`, which stores the result in
`As[1]`
* `get_bsxfun(op::Function)` and `get_bsxfun!(op::Function)` to instantiate a
broadcasting operation from a given per-element operation
* `bsx_getindex(A::AbstractArray, inds::Array...)`,
a pointwise broadcasting indexing operation as described in #2591,
and a matching `bsx_setindex!(A::AbstractArray, X::Array, inds::Array...)`

It should be straightforward to extend the functions to work with `StridedArray` arguments instead of just `Array`, but I'm uncertain of the proper interface to work with subarrays using manual linear indexing.

The code generates inner loop functions such as

function bsx_plus_inner!(shape::NTuple{2,Int}, arrays::NTuple{2,Array},
strides::Matrix{Int}, result::Array)
@assert (size(strides)==(2,2))
(n1,n2) = shape
if (n1==0); return; end
if (n2==0); return; end
(A1,A2) = arrays
k1 = 1
k2 = 1
(s1_1,s2_1,s1_2,s2_2) = strides
k = 1
for i2 = 1:n2
for i1 = 1:n1
begin
result[k] = +(A1[k1],A2[k2])
k += 1
end
k1 += s1_1
k2 += s2_1
end
k1 += s1_2
k2 += s2_2
end
end

The example above is for `bsxfun(+, A, B)` with two-dimensional result.
A linear index is updated using strides for each array, which should allow to accomodate
both broadcasting and subarrays.

On my machine, the timing example from #2991 gives

elapsed time: 0.806272502 seconds # subcol_forloop(a, b)
elapsed time: 3.1062779 seconds # subcol_bsxfun(a, b)
elapsed time: 0.936735113 seconds # .- from this pull request

i.e. the generated `.-` is only 16% slower than the handcoded example in #2991.

It's a somewhat complex piece of code (that generates code that generates code), and I think that it warrants some discussion. I could polish this forever, but I think that I need some reactions before moving forward. What do you guys think? Would something like this be useful? What parts of it?