Add operator simplification for unary negation and flatten nested AddedOperators

[albertomercurio]

Checklist

• Appropriate tests were added
• Any code changes were done in a way that does not break public API
• All documentation related to code changes were updated
• The new code follows the
  contributor guidelines, in particular the SciML Style Guide and
  COLPRAC.
• Any new documentation only uses public API

Additional context

I add some specific methods for unary negation that may significantly simplify the expression when applicable.

Moreover, I have added a function that flattens the input Tuple of an AddedOperator. This avoids to have nested AddedOperators when just calling AddedOperator(ops).

If this pr will be merged, could a new version be released?

[ChrisRackauckas]

this allocates though

[albertomercurio]

I agree, but this operation is usually performed once in the beginning when defining the operator. Then, the final operator is usually applied many times when solving an ODE or a linear problem.

This allows to have a more clean structure of the operator, which will be used repeatedly on a solver.

[ChrisRackauckas]

In the ODE it will allocate for every step which takes a new Jacobian.

[albertomercurio]

Could you make an explicit example? In my mind this is my workflow

some_f(a,u,p,t) = ...

# The allocations apply here, if the matrix is constant
A = ScalarOperator(rand(), update_function=some_f) * MatrixOperator(...) + ....

prob = ODEProblem(A, u0, tspan)
sol = solve(prob, Tsit5()) # It just updates the ScalarOperators. No allocations

[ChrisRackauckas]

that's not using the WOperator

[ChrisRackauckas]

You'd have to do a case where someone directly writes jac_prototype = MatrixOperator(J) I think, then it should allocate the extra terms instead of doing the lazy computation that is expected.

[albertomercurio]

Ok now I'm doing

f = ODEFunction(brusselator_2d_loop; jac_prototype = MatrixOperator(float.(jac_sparsity)))
prob_ode_brusselator_2d_sparse = ODEProblem(f, u0, (0.0, 11.5), p)

function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
    if newW === nothing || newW
        Pl = IncompleteLU.ilu(convert(AbstractMatrix, W), τ = 50.0)
    else
        Pl = Plprev
    end
    Pl, nothing
end

solve(prob_ode_brusselator_2d_sparse,
    KenCarp47(; linsolve = KrylovJL_GMRES(), precs = incompletelu,
    concrete_jac = true); save_everystep = false);

But still the same performances.

[ChrisRackauckas]

And memory usage and memory peak? The issue I'm worried about is that -A not being in-place will make there be a new temporary, and so it doubles the memory requirement on the system. I think this will oom larger cases.

[albertomercurio]

I'm actually trying

Base.:-(L::AbstractSciMLOperator) = (println("HELLO"); ScaledOperator(-true, L))

without the other method definitions, just to check that it is actually calling it. But this is not. So maybe we should try another case?

[albertomercurio]

Ok I have removed all the definitions that should allocate a new MatrixOperator, just to make it allocation-free.