Singular MM's initialization

[ChrisRackauckas]

Fixes #1039
Fixes #1042

@YingboMa can you add the singular mass-matrix initialization?

@kanav99 what should we do about du? Re-compute it on demand when we reinitialize?

P.S. need a doc update that describes the intializealgs

[coveralls]

Coverage increased (+3.1%) to 87.132% when pulling e79ebf6 on initialize into 1978905 on master.

[coveralls]

Coverage decreased (-19.0%) to 65.021% when pulling b4ab0e0 on initialize into 1978905 on master.

[codecov]

Codecov Report

Merging #1070 into master will decrease coverage by 1.05%.
The diff coverage is 81.21%.

@@            Coverage Diff             @@
##           master    #1070      +/-   ##
==========================================
- Coverage    77.2%   76.15%   -1.06%     
==========================================
  Files          95       95              
  Lines       31261    31448     +187     
==========================================
- Hits        24136    23950     -186     
- Misses       7125     7498     +373

Impacted Files	Coverage Δ
src/OrdinaryDiffEq.jl	100% <ø> (ø)	⬆️
src/integrators/integrator_utils.jl	89.23% <ø> (+1.17%)	⬆️
src/integrators/type.jl	100% <100%> (ø)	⬆️
src/caches/dae_caches.jl	86.2% <100%> (ø)	⬆️
src/misc_utils.jl	76.19% <100%> (+1.19%)	⬆️
src/integrators/integrator_interface.jl	61.29% <50%> (-4.73%)	⬇️
src/solve.jl	66.66% <72%> (-25.84%)	⬇️
src/initialize_dae.jl	81.76% <82.09%> (+35.93%)	⬆️
src/composite_solution.jl	19.04% <0%> (-28.58%)	⬇️
src/composite_algs.jl	68% <0%> (-28%)	⬇️
... and 29 more

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 1978905...e79ebf6. Read the comment docs.

[kanav99]

When all could the user ask for du? inside callbacks I guess and after each step! - any other case? In full initialization, du will always be set automatically. Otherwise in most BDFs, we have a linear formula of du we can have a separate function for each BDF. Note that du is never used internally and is not important for the algorithm to work.

[ChrisRackauckas]

I guess it only shows up in initialize, but with this method we need a du to reinitialize after a callback, so we somehow have to pick one.

[ChrisRackauckas]

I added initialization for singular mass matrices, but it assumes semi-explicit form.

[ChrisRackauckas]

Open an issue on semi-explicit form in DAE initialization

[ChrisRackauckas]

Open an issue on tolerances of DAE initialization

[ChrisRackauckas]

Talking with Yingbo, I realized that things that commonly use a non semiexplicit form are index 2, since they are relating differential variables to differential variables as a constraint. Seems like an issue to punt for now.

[kanav99]

Open an issue on semi-explicit form in DAE initialization

We don't have a problem form for Semi Explicit DAEs. Can I work on it next?

[dkarrasch]

This PR made the following problem fail, because it now calls lu on a Diagonal, which is not defined in LinearAlgebra:

using OrdinaryDiffEq, DiffEqDevTools, SparseArrays, LinearAlgebra

u0 = sin.(range(0, pi, length = 100))
tspan = (0.0, 1.0)

M = 0.5Diagonal(ones(100))
A = spdiagm(0 => ones(100))
update_func = (_A, u, p, t) -> _A.nzval .= t
f = DiffEqArrayOperator(A; update_func=update_func)

prob = ODEProblem(ODEFunction(f; mass_matrix=M), u0, tspan)
solve(prob, ImplicitEuler(linsolve=LinSolveFactorize(lu)))

gives

ERROR: MethodError: no method matching lu!(::Diagonal{Float64,Array{Float64,1}}, ::Val{true}; check=false)

Should there be specializations of issingular for structured matrices? Or maybe some keyword option, that allows to set issingular to true, but does an lu check by default?

[ChrisRackauckas]

Should there be specializations of issingular for structured matrices? Or maybe some keyword option, that allows to set issingular to true, but does an lu check by default?

Yes, we should specialize this on operator types and matrices.

[ChrisRackauckas]

Sorry if this is taking awhile, but I want to make sure we do this right. I'm taking this all the way to ArrayInterface to make issingular a thing: JuliaArrays/ArrayInterface.jl#42 and then we'll just be using that.