How to apply an in-place function to a slice of a vector?

New to Julia
, , , ,
1

I want to have an in-place function that takes a vector as an argument and applies another in-place function to a subarray (in general sense).
How do I do this?

function f!(du, u)
    println("\nf called")
    @show u
    du .= reverse(u)
    return nothing
end

"""
Applies f!() to the first 2 elements.
"""
function g!(du, u)
    f!(du[1:2], u[1:2])
    du[3] = 17
    return nothing
end

du = zeros(2)
u = [3, 11]

dug = zeros(3)
ug = [u; 0]

@show f!(du, u)
@show du

@show g!(dug, ug)
@show dug

gives

f called
u = [3, 11]
f!(du, u) = nothing
du = [11.0, 3.0]

f called
u = [3, 11]
g!(dug, ug) = nothing
dug = [0.0, 0.0, 17.0] # why not [11, 3, 17]?
2

In julia du[1:2] gives a copy of the slice, not a view (unfortunately). You can use the @view or @views macros, or the view function directly to construct views.

So your g! function should be

function g!(du, u)
    @views f(du[1:2], u[1:2])
    du[3] = 17
    return nothing
end
1 Like
3

That is, this creates a new object in memory, even without binding it to any name?

4

Bindings and allocations are totally different concepts, yes.

2 Likes
5

Does this behaviour affect only function calls? Do I understand correctly that slices passed by value?

In REPL, I can update the array with

u = zeros(3)
u[1:2] = ones(2)

and get [1, 1, 0] as expected.

6

No, we just made the (IMO misleading) choice to have u[1:2] = ones(2) to be special syntax for setindex!(u, ones(2), 1:2).

The u[1:2] there is never creating a slice, it’s just part of the setindex! syntax.

1 Like
7

I’m afraid I don’t get it completely.
I’m okay with the fact it’s some “special syntax”. But why does it work in REPL and not in a function (where it requires @view}? Where else will it not work?

8

the distinction isn’t repl vs function it’s before it after an equal sign.

6 Likes
9

When Julia sees arr[indices...] = v, it translates it to setindex!(arr, v, indices...) at an early stage. Docs: setindex!.

5 Likes
10

Incidentally, I believe reverse(u) will allocate a intermediare temporary vector. A lazy non-allocated version can be had by using Iterators.reverse or reverse the indices of a view, ie

du .= @view u[end:-1:begin]
4 Likes
11

Thanks. Reverse is there only to provide a MWE.

12

Sure, but since the question concerns in-place operations and allocations, and since I expect almost everything in an MWE to be relevant, I assumed this was too.

1 Like