Refactoring iterative interface solvers

First steps to make the addition of custom solvers easier and more modular. Adds an abstract class describing a solver whose interface is a solve() method taking the system as input. Currently, only Jacobi is implemented as an illustration. If Jacobi isn't used, we default to PETSc solvers as previously.

Open design questions:

• Who should own the residual history ? Currently the formulation does, but I think it would make more sense to have the solver own it. Maybe move this to the common Solver abstract class ? Alternative would be to have the solver "feed" the formulation by telling it the residual at each iteration. Owning it could help implement things like "pick restarts depending on previous residuals" etc.
• Should tolerance / max iter be an argument of Formulation::solve() (as currently), or be a state of the solver object ? In the long run I'd probably refactor Formulation::solve() to take a Solver object that knows these things, as well a solver-specific parameters (e.g. GMRES restart values)
• How to handle Krylov solvers from PETSc (KSPxx) ? For now they are simply matched by strings. Should we do a simple class "PETSc external solver" or something ?
• What happens to the _IA variable ? Should a solver set it ? Should ownership change from the formulation to the solver ? What does it do exactly ? I see it changes the output of a matrix-vector multiplication by removing one time X, so it computes (A-I)*x ? I'm slightly confused because we usually say we solve (I-A)x=b and the sign convention is not fully familiar to me.

Also, what was the reason why the function MatVectProduct was static ? Since it's not part of a class this seems like C-style code, and it generates a warning I don't fully understand, so I removed it here. Should it become a static member of the Formulation class ?