Previous attempts to integrate SciMLOperators into the ecosystem broke everything. This was because we coupled together two separate tasks:

1. supporting SciMLOperators
2. deprecating AbstractDiffEqOperator, DiffEqLinearOperator, AbstractDiffEqCompositeOperator, and the concrete types.

So we were never able to get anywhere because errors (in downstream packages) from the deprecations in step 2 were too many to handle. This time we shall do the integration step-by-step.

• Define all operators and traits in SciMLOperators
• Have SciMLBase reexport SciMLOperators. Do NOT deprecate diffeqoperator family
• Ditto for DiffEqBase.jl
• Add support for SciMLOperators in downstream packages.
• Replaces subtypes of *DiffEqOperators to AbstractSciMLOperators. Do NOT remove support for DiffEqOperators. Order below:
• LinearSolve.jl - make preconditioners AbstractSciMLOperators
• Sundials.jl
• DiffEqOperators.jl
• SparseDiffTools.jl - make AD operators AbstractSciMLOperators
• DiffEqSensitivity.jl - use SparseDiffTools AD operators in SteadyStateAdjoint. Use adjoint(A) if islinear(A) = true, has_adjoint(A) = true.
• OrdinaryDiffEq.jl - form WOperator just by composing AbstractSciMLOperator: W = 1/gamma * M - J or the opp
• StochasticDiffEq.jl
• Once SciMLOperators are supported everywhere, remove DiffEqOperators starting from downstream packages to upstream packages.
• Finally deprecate AbstractDiffEqOperator family in SciMLBase.jl
• remove DiffEqOperators from docs

EDITS -

TODOs

SciMLBase.jl

deprecate all diffeqoperator subtypes

LinearSolve.jl

• update preconditioner docs
• deprecations:

const ComposePreconditioner = SciMLOperators.ComposedOperator
@deprecate ComposePreconditioner SciMLOperators.ComposedOperator

const InvPreconditioner = SciMLOperators.InvertedOperator
@deprecate InvPreconditioner SciMLOperators.InvertedOperator

• in defaultalg selection, as well as in init_cacheval, OperatorAssumptions{nothing} behaves exactly like OperatorAssumptions{true}. So let's remove the option to put in nothing.
• in src/common.jl, set_A also modify cache.OperatorAssumptions.
• in src/default.jl, remove method defaultalg(A::Diagonal, b, ::OperatorAssumptions{false}) because diagonal matrices are always square
• In src/factorization.jl, DiagonalFactorization should ideally preinvert the D.diag in init_cacheval, and then mul! in SciMLBase.solve. Create a separate PR
• replace IterativeSolvers: Identity in OrdinaryDiffEq with IdentityOperator.

DiffEqBase.jl

In test/basic_operator_tests.jl, test/affine_operator_tests.jl, there are tests for DiffEqOperators. Let's leave that as it is till DiffEqOps are remove from downstream packages.

SparseDiffTools.jl

• Add OrdinaryDiffEq.jl - InterfaceII to downstream tests. Same for SciMLSensitivity.jl when we get done with that.
• Add JacVec ODE tests
• autodiff shouldn't be boolean logic but multi-valued logic through types choices, like AutoZygote(). change the autodiff choices to be non-boolean logic and Zygote an extension package in order to release.
• release v2 Release v2 JuliaDiff/SparseDiffTools.jl#230

# JacVec OrdinaryDiffEq itnegration test

function lorenz(du, u, p, t)
    du[1] = 10.0(u[2] - u[1])
    du[2] = u[1] * (28.0 - u[3]) - u[2]
    du[3] = u[1] * u[2] - (8 / 3) * u[3]
end
u0 = [1.0; 0.0; 0.0]
tspan = (0.0, 100.0)

ff1 = ODEFunction(lorenz, jac_prototype = JacVec(lorenz, u0))
ff2 = ODEFunction(lorenz, jac_prototype = JacVec(lorenz, u0, autodiff=false))

for ff in [ff1, ff2]
    prob = ODEProblem(ff, u0, tspan)
    @test solve(prob, TRBDF2()).retcode == :Success
    @test solve(prob, TRBDF2(linsolve = KrylovJL_GMRES())).retcode == :Success
    @test solve(prob, Exprb32()).retcode == :Success
    @test solve(prob, Rosenbrock23()).retcode == :Success
    @test solve(prob, Rosenbrock23(linsolve = KrylovJL_GMRES())).retcode == :Success
end