Should an in-place FunctionOperator be totally allocation free?

General Usage
1

The code

    fop = FunctionOperator((v, u, p, t) -> mul!(v, K_ff, u), F_f, zeros(length(F_f)))
    prob = LinearProblem(fop, F_f)
    @time sol = solve(prob, ALG(), Pl=PRECOND, verbose=verbose)

reports allocations:

  7.333552 seconds (966 allocations: 25.764 MiB)

More than if I use MatrixOperator.

    mop = MatrixOperator(K_ff)
    prob = LinearProblem(mop, F_f)
    @time sol = solve(prob, ALG(), Pl=PRECOND, verbose=verbose)

yields

  7.140111 seconds (42 allocations: 25.693 MiB)

I am guessing the allocations must be some small change, given that the amount of allocated memory is nearly the same? But: is this as expected or am I doing something wrong?

2

Could you please post some complete working code, so others can experiment?

3

To me, this smells like the usual issue of boxing in closures.

Try e.g.

fop = let K  = K_ff, F = F_f
FunctionOperator((v, u, p, t) -> mul!(v, K, u), F, zeros(length(F)))
end

But ideally, don’t use closures – it sucks to audit them for boxing. Or are there good linters that reliably warn about closures that box? Then enable that for your project.

Instead define a nice type

struct Wrapper{KT, FT} <: Function
K::KT
F::FT
end
(w::Wrapper)(v, u, p, t) = (mul!(v, w.K, u), w.F, zeros(length(w.F)))
fop = FunctionOperator(Wrapper(K_ff, F_f))

In my view, capture-by-binding is a total language misfeature in julia; ideally we’d have capture-by-value.

Luckily you can have capture-by-value semantics in a well-defined subset of julia – you just need to stick to only capturing SSA values. Ideally with help of a linter that complains whenever you capture something that is not SSA.

1 Like