Broadcasting question

New to Julia
1

I’m trying to calculate the following:
image

where v and \hat{n} vectors of length 3. But I need to do this to a bunch of vectors stored in a Vector{Vector{SVector{3,Float64}}}

I ended up doing it like this:

dots = [[velocities[j][i]' * nHat[j][i] * nHat[j][i] for i =1:length(config.atomicBasis[j])] for j = 1:config.order]
velocities = velocities - 2 * dots

but wondering if there is a better way.

Thanks in advance

2

You’re allocating a bunch of temporary arrays, e.g. dots = [...] is allocating lots of little arrays, 2 * dots allocates new arrays, then velocities - 2 * dots allocates more arrays.

Assuming you want to update velocities in-place, it would be a lot clearer (and faster, with no temporary arrays) to just write a loop:

for j = 1:config.order, i = 1:length(config.atomicBasis[j])
   @inbounds velocities[j][i] -= (2 * (velocities[j][i]' * nHat[j][i])) * nHat[j][i]
end

(Hopefully the compiler can hoist the velocities[j] and nHat[j] indexing out of the innermost loop?)

Loops (in local scope) are fast in Julia. Don’t twist yourself into a pretzel trying to use “vectorized” operations.

3 Likes
3

You could take the advice of stevengj (which is what you actually should do), or you could go julia-nuts:

using StaticArrays
using LinearAlgebra

new_vec(v::SVector, n̂::SVector) = v - 2*(v⋅n̂)*n̂
new_vec(v::Vector{<:SVector}, n̂::Vector{<:SVector}) = new_vec.(v, n̂)
new_vec(v::Vector{<:Vector{<:SVector}}, n̂::Vector{<:Vector{<:SVector}}) = new_vec.(v, n̂)
new_vec!(v::Vector{<:SVector}, n̂::Vector{<:SVector}) = v .= new_vec.(v, n̂)
new_vec!(v::Vector{<:Vector{<:SVector}}, n̂::Vector{<:Vector{<:SVector}}) = v .= new_vec!.(v, n̂)

vs = [rand(SVector{3, Float64}, 4) for _ ∈ 1:2]
n̂s = [rand(SVector{3, Float64}, 4) for _ ∈ 1:2]

vs = new_vec!(vs, n̂s)
5

I’ve noticed this theme emerge for me: I’ve been coding in python too long and always think I should write complicated code in the fewest number of lines possible. Would it be fair to say that the more julian thing do is often to just write out the dang loop?

6

Yes, Python trains you that your code is slow, and libraries are fast, so you should contort yourself to make all loops occur in libraries. To some extent, this is good advice in any language — re-using highly optimized, robust libraries is good!

But it has two limitations. First, you sometimes need to write really convoluted, opaque code to express your problem in terms of the library API. Second, you often pay a performance price for using a sequence of generic library functions instead of code that is specialized for your problem. (For example, you might need several “vectorized” calls instead of a single loop, which slows things down a lot. Julia’s “dot call” loop fusion can help with this, but not so much with nested data structures like yours.)

So, it’s a judgement call, but in Julia “writing out the dang loop” is a good alternative to writing convoluted “vectorized” code.

3 Likes