Operator mass matrices: `IdentityOperator`/`ScalarOperator` fail at solve entry — fix in place or push toward native operator support?

[singhharsh1708]

IdentityOperator and ScalarOperator cannot be used as f.mass_matrix right now because ArrayInterface.issingular isn't defined for them, so the DAE check at lib/OrdinaryDiffEqCore/src/solve.jl:233 throws before the solver even starts.

Repro against master (dev'd sublibs):

using OrdinaryDiffEqSDIRK, SciMLOperators, LinearAlgebra

f(du, u, p, t) = (du .= -u; nothing)
for mm in (IdentityOperator(3), ScalarOperator(2.0))
    prob = ODEProblem(ODEFunction(f, mass_matrix = mm), ones(3), (0.0, 1.0))
    try
        solve(prob, TRBDF2(); reltol = 1e-6, abstol = 1e-8)
    catch e
        println(typeof(mm).name.name, " => ", first(sprint(showerror, e), 120))
    end
 end

IdentityOperator => MethodError: no method matching issingular(::IdentityOperator)
ScalarOperator   => MethodError: no method matching issingular(::ScalarOperator{Float64, ...})

MatrixOperator and DiagonalOperator (which is a MatrixOperator under the hood) work fine because build_J_W has a dedicated MatrixOperator branch.

The third TODO on top of build_J_W (lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl:1100-1102) already anticipates this:

# TODO - make J, W AbstractSciMLOperators (lazily defined with scimlops functionality)
# TODO - if jvp given, make it SciMLOperators.FunctionOperator
# TODO - make mass matrix a SciMLOperator so it can be updated with time. Default to IdentityOperator

So two directions this could go, and I would rather hear which one you want before I open PRs:

1. Point fix: add ArrayInterface.issingular methods for IdentityOperator / ScalarOperator (and any other operator that could legitimately appear as a mass matrix) upstream in SciMLOperators or ArrayInterface. Keeps concrete matrices as the internal representation.

2. Native operator direction (matches the TODO): promote f.mass_matrix to an AbstractSciMLOperator at some boundary and collapse the mass_matrix === I / isa UniformScaling / <matrix> branches to single operator dispatch. Bigger change, closes the whole TODO.

If (2): should the promotion happen at ODEFunction construction in SciMLBase, or inside each subpackage's alg_cache?

[ChrisRackauckas]

Point fix: add ArrayInterface.issingular methods for IdentityOperator / ScalarOperator (and any other operator that could legitimately appear as a mass matrix) upstream in SciMLOperators or ArrayInterface. Keeps concrete matrices as the internal representation.