A .= circshift(b, shfts) and circshift!(a, b, shfts)

Hi all,

I’m learning Julia by reimplementing some of my numerical algorithms.

As far as I understand from the manual and language design, the following two operations should be equivalent:

a .= circshift(b, shifts)
    0.422910 seconds (13 allocations: 762.940 MiB)

and

circshift!(a, b, shifts)
    0.164945 seconds (11 allocations: 432 bytes)

However it looks that they are executed differently and memcpy takes place in .= variant. Am I missing something, or is it a bug\language inconsistency?

Thanks.

1 Like

The first version results in an allocated array which entries are then in-place assigned to a. Therefore there is an allocated array somewhere within the circshift function which consumes memory.

2 Likes

Sort of rephrasing what @roflmaostc pointed out:

circshift doesn’t know that there will be an in-place assignment afterwards (i.e. it can not be “fused”).

circshift!, however, does know that there is already preallocated memory which can be used.

2 Likes

You can use ShiftedArrays.circshift for a lazy, non-allocating version of circshift (see docs):

julia> using ShiftedArrays, BenchmarkTools

julia> a, b = rand(100), rand(100);

julia> @btime $a .= ShiftedArrays.circshift($b, 1)
  128.619 ns (0 allocations: 0 bytes)
2 Likes

Thank for replies. I understand what is going on. Standard circshift has to create a buffer. But here https://docs.julialang.org/en/v1/manual/functions/#man-vectorized it is stated that X .= ... allows to preallocate output. Which is not followed for circshift. Does output preallocation works only for trivial math?

do you know if these view abstractions will work efficiently with ffts?

In ffts the output will be indeed a ShiftedArrays view and hence that in-place update should work.
However, fft itself allocates intermediate arrays.
There is also ffts!

See my previous reply. But in this case, what is the whole idea of X.=... syntax? It seems to be trivial that if two implementation are given, e.g:
f(x)
and
f!(out, x)
then z.=f(x) can be evaluated as f!(z, x)?

ffts!” ? Did you mean inplace fft or efficient fft with shift? I couldn’t find it in AbstractFFTs.jl

The dot essentially means “element-wise”. So z .= ... means write ... element-wise into x (think of a for loop). So, a typical use case is something like z .= f.(x) where the operations “apply f element-wise to x” and “assign the result of f.(x) element-wise to z” get fused into a single “for loop” so to avoid a temporary allocation for f.(x).

If you don’t do element-wise applications, like you do with circshift, the dot syntax doesn’t help much (as you’ve seen in your first case). That’s why dedicated in-place functions like circshift! exist.

Hope that clears things up a bit?

3 Likes

It might be instructive to compare and understand the difference between

x .= a .+ b # x .= f.(x)

and

x .= a + b # x .= f(x) (no broadcasting of `f`)
1 Like

Ok, so x .= is intended to work only for trivial element-wise operations only? I.e. in your example, if I had a manually vectorized function f(x) then
z .= f(x)
would copy data, but autovectorized
z .= f.(x)
would do everything inplace?

Oh sorry, I thought by “ffts” you are referring to a package I maintain :smiley:

Look, here at FourierTools.jl

The thing is that, in contrast to MATLAB (as far as I remember), you wouldn’t want to define a vectorized version of f, say f(y::Matrix{Float64}), but only a scalar variant f(x::Float64) to then simply use f.(y).

Example:

# Let's say we want to compute the `exp(sin(x))` of all elements of a matrix
f(x::Float64) = exp(sin(x))

M = rand(10,10) # our matrix
R = zeros(10,10) # preallocated matrix to store the result

R .= f.(M)


1 Like

Well, whatever you mean by “trivial” :grinning_face_with_smiling_eyes:

But yes, there is no generic magical way in Julia (or any other language?) to autogenerate in-place versions of “non-trivial” matrix functions.

1 Like

The thing is, in my code I “vectorize” over matrices. But thanks for your help! I think I see the logic behind the syntax. Only comment I have, is that docs ulitmately state that .= allows to preallocate input, but this mechanism is actually very limited, as we just discussed.

If I may ask another question. Is there any standard macro that can translate z .= f(x) to f!(z, x) ? It would be very helpful for algorithm readability.

Very cool, thanks for the reference. Are you planning to add inplace ft! ? That’s exactly what I would need as well :))

Well, either use dot-broadcasting then or define in-place versions of your functions. But maybe it isn’t even necessary to “vectorize over matrices”? Perhaps a simple loop will do as well? (I don’t know. Depends on context of course.)

.= doesn’t preallocate anything. It writes into preallocated memory if you will. Where is this mentioned in the docs? Can you point out the paragraph / line?

Well, as we’ve tried to explain above, these two things do very different operations, so why would you want a macro to translate one to the other? Also, not for every f there is an in-place pendant f!.

ft! is not possible since fft! can’t act in place on a ShiftedArrays view. Hence FFTW.jl needs to collect the array.
The performance drawback is usually not that critical but if you run out of memory that’s of course an issue.