Improve speed and reduce allocations in the ManifoldUpdate callback

[nathanaelbosch]

Currently, the manifoldupdate function is implemented very inefficiently. By using the cache matrices already present in integ.cache, it should be possible to implement this in a much more efficient manner.

[nathanaelbosch]

In large parts this is fixed in #134. But there is still room for more improvements! Allocations could be further reduced by

• requiring that the residual function is in-place
• assuming that the residual has the same shape as the ode solution u
• generally re-using already existing cache matrices; this is easies in conjunction with the points above

The issue is that, when we assume that the residual has the same dimensionality as the ODE, it often has only one entry and the rest are zeros. This leads to low-rank Jacobians, and to low-rank measurement covariances, and then we run into singularity errors. So, the main thing to be figured out here next is:
How can we (efficiently) implement the Kalman update such that it can handle singular measurment covariances?

Some example code to work on this issue:

using ProbNumDiffEq

function harmonic_oscillator(du, u, p, t)
    du[1] = u[2]
    du[2] = -u[1]
end
prob = ODEProblem(harmonic_oscillator, ones(2), (0.0, 10.0))

E_old(u) = [dot(u, u) - 2]
E_new(u) = [dot(u, u) - 2; 0]

solve(prob, EK1(), callback=ManifoldUpdate(E_old)); # runs
solve(prob, EK1(), callback=ManifoldUpdate(E_new)); # throws a SingularException

The desired outcome here is that E_new runs, while E_old fails due to incorrect shapes. Then, in a separate next step we can look into in-pace residual functions and re-using our cache matrices to remove the remaining allocations.