Why does `./=` allocate?

General Usage
1

I have this bit of code. I managed to avoid all allocations.

function compute!(csmatout, XYZ, tangents, feid, qpid)
        center = (0.0, 0.0, 0.0)
        xyz = (XYZ[1], XYZ[2], XYZ[3])
        csmatout[:, 1] .= xyz .- center
        csmatout[3, 1] = 0.0
        # csmatout[:, 1] ./= norm(@view(csmatout[:, 1]))
        nm = norm(@view(csmatout[:, 1]))
        for k in 1:3; csmatout[k, 1] /= nm; end
        csmatout[:, 3] .= (0.0, 0.0, 1.0)
        cross3!(@view(csmatout[:, 2]), @view(csmatout[:, 3]), @view(csmatout[:, 1]))
        # csmatout[:, 2] ./=  norm(@view(csmatout[:, 2]))
        nm = norm(@view(csmatout[:, 2]))
        for k in 1:3; csmatout[k, 2] /= nm; end
        return csmatout
end

But, I had to replace the lines that are commented out. They would allocate. In other words, this allocates.

function compute!(csmatout, XYZ, tangents, feid, qpid)
        center = (0.0, 0.0, 0.0)
        xyz = (XYZ[1], XYZ[2], XYZ[3])
        csmatout[:, 1] .= xyz .- center
        csmatout[3, 1] = 0.0
        csmatout[:, 1] ./= norm(@view(csmatout[:, 1]))
            # nm = norm(@view(csmatout[:, 1]))
            # for k in 1:3; csmatout[k, 1] /= nm; end
        csmatout[:, 3] .= (0.0, 0.0, 1.0)
        cross3!(@view(csmatout[:, 2]), @view(csmatout[:, 3]), @view(csmatout[:, 1]))
        csmatout[:, 2] ./=  norm(@view(csmatout[:, 2]))
            # nm = norm(@view(csmatout[:, 2]))
            # for k in 1:3; csmatout[k, 2] /= nm; end
        return csmatout 
end

I thought that the compiler would be able to divide the column without allocating. What is it that I’m misunderstanding here?

Thanks!

1 Like
2

csmatout[:, 1] ./= nm is equivalent to

csmatout[:, 1] .= csmatout[:, 1] ./ nm

so the csmatout[:, 1] on the right-hand side allocates. This can be fixed by putting @views before the lines, but it will be a lot simpler to just to put @views on the whole function and get rid of all the @view calls.

i.e. something like:

@views function compute!(csmatout, XYZ, tangents, feid, qpid)
        center = (0.0, 0.0, 0.0)
        xyz = (XYZ[1], XYZ[2], XYZ[3])
        csmatout[:, 1] .= xyz .- center
        csmatout[3, 1] = 0.0
        csmatout[:, 1] ./= norm(csmatout[:, 1])
        csmatout[:, 3] .= (0.0, 0.0, 1.0)
        cross3!(csmatout[:, 2], csmatout[:, 3], csmatout[:, 1])
        csmatout[:, 2] ./=  norm(csmatout[:, 2])
        return csmatout
end

Isn’t it more readable without all of those @view calls?

PS. If you are working with lots of 3-component vectors representing XYZ coordinates (so that you know their length at compile time), I would strongly consider arrays of SVector from StaticArrays.jl instead.

12 Likes
3

Ah, I forgot I should use @view for the lhs too.
But @views sure is so much nicer!
Ta!

4

This always gets me and I guess I’m not alone. I also forgot why we can’t make

ary[n:m] ./= 3

lowers to equivalent to having a view.

Actually, now I remember, THIS REALLY ALWAYS gets me:

ary1[n:m] .= @view ary2[n:m]

is enough but not for ary1[n:m] .*=

1 Like