GMRES tie-ins for linsolve support

[ChrisRackauckas]

@isaacsas has mentioned that there are some extra tie-ins we will want to have for full GMRES support. Specifically, we need a way to specify how the initial condition of the GMRES should be done. Right now it's hardcoded to always start at 0: https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src/linear_nonlinear.jl#L60 . For our NLNewton we should make sure that this vector is instead something that makes sense for an initial GMRES guess, either from the last stages or zero'd. The Rosenbrock methods need to handle this as well. This was explicitly put in because in some cases it was the direct output of similar which of course is unpredictable.

Another thing that comes up is whether the JacVecOperator should be updating along with the J*v calculations. We probably just need to better specify when to update coefficients:

https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/src/jacvec_operators.jl#L93-L130

to have that done automatically. That might be done in the GMRES via a tie-in to the preconditioner?

But the last thing is that in this formulation we don't cache the Krylov subspace. Is it a good idea to do so? Someone should dig into the Sundials code and find out. If it is a good idea, is there another Julia package we should use here? Is there an easy modification to IterativeSolvers? Should we roll our own?

@YingboMa

[YingboMa]

Another thing that we need to tone for GMRES is that it can give you exact Newton iteration, so the convergence rate is quadratic.

[isaacsas]

The comment by @YingboMa is what I was trying to get at on Slack. The Krylov solver approaches that are matrix free would benefit from using an inexact Newton method which updates W with the current solution estimate during each step within a nlsolve call. I suspect this is the cause of the very poor GMRES performance on the BCR network.

[ChrisRackauckas]

That's what I meant by:

Another thing that comes up is whether the JacVecOperator should be updating along with the J*v calculations. We probably just need to better specify when to update coefficients:

https://github.com/JuliaDiffEq/DiffEqOperators.jl/blob/master/src/jacvec_operators.jl#L93-L130

to have that done automatically. That might be done in the GMRES via a tie-in to the preconditioner?

We just need to stick an update_coefficients in there somewhere. Can we just hack the preconditioner to do update coefficients before J*v? Or do we need to roll our own for other reasons?

[isaacsas]

Sorry, I missed that. Yes, an update_coefficients is what is needed. I guess we would want to have this done in a way that a user-supplied preconditioner to GMRES still works ok too. Would it be problematic to just have a flag or trait associated with the linear solver that causes nlsolve! to just call this?

More generally, maybe one might want to use different stopping criteria for an inexact Newton method vs modified Newton? So a longer-term goal might be to get an inexact Newton solver (once swapping nonlinear solvers is better supported). For now just getting an update_coefficients call into the Newton solver might still give a big performance boost.

[ChrisRackauckas]

I guess we would want to have this done in a way that a user-supplied preconditioner to GMRES still works ok too.

By the hack, I mean make a preconditioner which does update_coefficients and acts like I, but then handle a user passed in preconditioner by adding it to the anonymous function we are making so that it does update_coefficients and then acts how the user expects. That's why it's a bit hacky, but it would work :).

Would it be problematic to just have a flag or trait associated with the linear solver that causes nlsolve! to just call this?

Or we could always do it in nlsolve! whenever it's a DiffEqOperator. That's probably a better solution.

More generally, maybe one might want to use different stopping criteria for an inexact Newton method vs modified Newton? So a longer-term goal might be to get an inexact Newton solver (once swapping nonlinear solvers is better supported).

We need a lot more tie-ins than we currently have. I think making a more solid interface for these things might be required down the line.

[isaacsas]

Maybe JacVecOperator should be an AbstractMatrixFreeOperator, and that is what the check is over? Then explicit matrix-based DiffEqOperators will still use the modified Newton method.

[ChrisRackauckas]

Complete.