Dual linear solvers interface

[jagot]

In relation to #90, and through discussions on Gitter, there is an idea of creating a dual interface to the solvers for linear problems:

1. For standard methods such as implicit Euler and the trapezoidal rule (Crank–Nicolson), for each DiffEqOperator at least a A_mul_B! should be provided, and optionally a \.

   • For certain types, \ can be inferred, e.g. for Tridiagonal, lufact! could be used. If the operator is not time-dependent (except for an overall scalar coefficient), the lufact! can be performed ahead of time.
   • For generic matrix-free operators (as well as sparse matrices, etc), a solver such as GMRES could be inferred instead.

2. Alternatively, the user might specify an optimized expmv! implementation.

The solution manner should be inferred from the operator, i.e. has_expmv, and it should be possible to mix-and-match, i.e. for u' = F(u,t)u with F = A + B, the user should be able to provide a expmv! for B and use a Crank–Nicolson propagator for A:

expmv!(tmp1, B, dt/2, u)
expmv!(tmp2, A, dt, tmp1)
expmv!(v, B, dt/2, tmp2)

where expmv!(tmp2, A, dt, tmp1) resolves through dispatch to a generic exponentiator (probably the user can set a flag in the DiffEqOperator that :CrankNicolson/:ImplicitEuler is to be used).

---

Side note: I am unsure whether the time-step should be absorbed into the operator or not; for Krylov approximations to the exponential, one usually computes the subspace without the time-step, and applies it in the subspace exponential (notation the same as in http://epubs.siam.org/doi/10.1137/0729014):

exp(τA)v ≈ βVₘₘ_m exp(τH_m)e₁

where H_m is the approximation of A in the mth Krylov subspace spanned by the columns of V_m.

[ChrisRackauckas]

Side note: I am unsure whether the time-step should be absorbed into the operator or not; for Krylov approximations to the exponential, one usually computes the subspace without the time-step, and applies it in the subspace exponential

That can be the choice of the solver. If the solver does a direct t*A it absorbs. But it can write the algorithm other ways of course.

[ChrisRackauckas]

has_expmv, has_expm, has_mul, has_ldiv should be traits added to the interface that the methods can query for a given operator to throw errors or try something different.

where expmv!(tmp2, A, dt, tmp1) resolves through dispatch to a generic exponentiator (probably the user can set a flag in the DiffEqOperator that :CrankNicolson/:ImplicitEuler is to be used).

Hmm, I'm not so sure about this. I'll have to see what an implementation looks like.

[jagot]

Good reference on different exponential propagators (applied to QM, but written by a mathematician):

@Book{lubich2008from,
 author = {Lubich, Christian},
 title = {From quantum to classical molecular dynamics : reduced models and numerical analysis},
 publisher = {European Mathematical Society},
 year = {2008},
 address = {Zürich, Switzerland},
 isbn = {978-3-03719-067-8}
}

[ChrisRackauckas]

We now have a good linear solver interface.