Help with coding up some matrix / broadcast function

General Usage
1

Say I’ve coded up a type called Field, and defined broadcast over these objects so that I can multiply them, i.e.

f1 = Field()
f2 = Field()

f1 .* f2 .* f1 # this works correctly

Now I would also like to be able to work with Vectors of Fields, in particular I want matrix multiplication to “use the broadcasted * operator”, and to do so efficiently. By that I mean,

x = [f1,f2]'
y = [f1,f2]

# I want to be able to write this:
x * y + f1
# and have this be effectively evaluated:
x[1] .* y[1] .+ x[2] .* y[2] .+ f1

In particular, I want to make sure there’s one single broadcast over everything, so that even the x * y + f1 form does not allocate any temporary arrays.

Can anyone think of a good way to do this?

My best idea so far involves writing a custom @generated broadcast function, and then always writing @. x * y + f1. From inside the generated broadcast function I then have access to both the expression x * y + f1, and the types of the arguments, x, y and f1. Hence in theory I have everything I need to transform the call into what is needed. However, I’m not the most keen on what it might take to code up the logic for this transformation so it will work generally. Is there perhaps an easier way?

Thanks for any advice.

2

I’m not sure how to implement the approach you’re envisioning. You can always define your own method of the * function for arrays with elements of type Field. For example (in Julia 0.6) you could :

struct Field
    a::Real
    b::Real
end

import Base.*
function *(A::Array{Field,2},B::Array{Field,2})
   #Implement your multiplication here
end

Then inside your new * function you can define multiplication of arrays with element type Field using broadcasting in the way you like. Is there a reason why that approach wouldn’t work for you?

3

Thanks, yea, the only reason that approach wouldn’t work is I’m trying to avoid allocating any temporary arrays and just use broadcasting across the entire expression. If I define * in the way you suggest, then x*y can certainly be done with a single broadcasted expression, but then the result is stored as a temporary array before being added to f1, rather than the whole thing being evaluated in one loop.

4

Why can’t you just write x .* y .+ f1, or @. x * y + f1? This will work in 0.6 without requiring you to define any new broadcast methods.

Fusion of x * y + f1 is not going to happen (at least, not in the near future), because it is a much harder problem — unlike the x .* y .+ f1 fusion, which can happen at parse time (really, “lowering” time), fusing x * y + f1 would have to happen much later, at compile time (after all of the types are identified etc.). This is the whole point of the dot syntax in Julia.

5

(Note I edited my original post where I had an x instead of a y that might have caused some confusion, sorry about that)

The reason is just that that’s a different operation than what I’m after. To try and explain a bit better, in my example of x*y+f1, x is a row vector of Fields, y is a column vector of Fields, and f1 is a Field, so I want the x*y part to represent a dot product like it normally does for row/column vector multiplication, not a broadcast. That is, I want that to effectively evaluate x[1] * y[1] + x[2] * y[2] + f1 but I want the whole thing to be broadcasted.

If I did x.*y.+f1 the answer is going to be a Matrix of Fields, so not the operation I’m after.