Solvers with KrylovJL fail

[DanielVandH]

Should the code below be expected to work?

using OrdinaryDiffEq, LinearSolve, Krylov, NonlinearSolve
function lorenz!(du,u,p,t)
    @show t
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan) |> NonlinearProblem
alg = NewtonRaphson(linsolve=KrylovJL_GMRES())
sol = solve(prob, alg)

ERROR: MethodError: no method matching forwarddiff_color_jacobian!(::Matrix{Float64}, ::NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, ::Vector{Float64}, ::Nothing)

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

  [6] #solve_call#22
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:509 [inlined]
  [7] solve_call
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:479 [inlined]
  [8] #solve_up#29    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:940 [inlined]  [9] solve_up    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:905 [inlined]
 [10] #solve#28
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:899 [inlined] [11] solve(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, args::NewtonRaphson{0, true, Val{:forward}, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing})    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:889
 [12] top-level scope
    @ Untitled-1:12

If instead autodiff = false:

alg = NewtonRaphson(linsolve=KrylovJL_GMRES(), autodiff = false)
sol = solve(prob, alg)

ERROR: MethodError: no method matching default_relstep(::Nothing, ::Type{Float64})

Some of the types have been truncated in the stacktrace for improved reading. To emit complete information
in the stack trace, evaluate `TruncatedStacktraces.VERBOSE[] = true` and re-run the code.

Closest candidates are:
  default_relstep(::Type, ::Any) at C:\Users\User\.julia\packages\FiniteDiff\e7xnl\src\epsilons.jl:25
  default_relstep(::Val{fdtype}, ::Type{T}) where {fdtype, T<:Number} at C:\Users\User\.julia\packages\FiniteDiff\e7xnl\src\epsilons.jl:26
Stacktrace:
  [1] finite_difference_jacobian!(J::Matrix{Float64}, f::NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, x::Vector{Float64}, fdtype::Nothing, returntype::Vector{Float64}, f_in::Nothing) (repeats 2 times)
    @ FiniteDiff C:\Users\User\.julia\packages\FiniteDiff\e7xnl\src\jacobians.jl:298
  [2] jacobian_finitediff_forward!(J::Matrix{Float64}, f::NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, x::Vector{Float64}, jac_config::Nothing, forwardcache::Vector{Float64}, cache::NonlinearSolve.NewtonRaphsonCache{true, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, NewtonRaphson{0, false, Val{:forward}, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, GmresSolver{Float64, Float64, Vector{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, true}, Matrix{Float64}, Nothing})
    @ NonlinearSolve C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\jacobian.jl:19
  [3] jacobian!(J::Matrix{Float64}, cache::NonlinearSolve.NewtonRaphsonCache{true, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, NewtonRaphson{0, false, Val{:forward}, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, GmresSolver{Float64, Float64, Vector{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, true}, Matrix{Float64}, Nothing})
    @ NonlinearSolve C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\jacobian.jl:45
  [4] perform_step!(cache::NonlinearSolve.NewtonRaphsonCache{true, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, NewtonRaphson{0, false, Val{:forward}, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing}, Vector{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, typeof(NonlinearSolve.DEFAULT_NORM), Float64, NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, NonlinearSolve.JacobianWrapper{NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, SciMLBase.NullParameters}, LinearSolve.LinearCache{Matrix{Float64}, Vector{Float64}, Vector{Float64}, SciMLBase.NullParameters, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, GmresSolver{Float64, Float64, Vector{Float64}}, LinearSolve.InvPreconditioner{LinearAlgebra.Diagonal{Float64, Vector{Float64}}}, LinearAlgebra.Diagonal{Float64, Vector{Float64}}, Float64, true}, Matrix{Float64}, Nothing})
    @ NonlinearSolve C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\raphson.jl:160
  [5] solve!
    @ C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\raphson.jl:188 [inlined]
  [6] #__solve#1
    @ C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\NonlinearSolve.jl:31 [inlined]
  [7] __solve
    @ C:\Users\User\.julia\packages\NonlinearSolve\Da5qG\src\NonlinearSolve.jl:27 [inlined]
  [8] #solve_call#22
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:509 [inlined]
  [9] solve_call
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:479 [inlined]
 [10] #solve_up#29
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:940 [inlined]
 [11] solve_up
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:905 [inlined]
 [12] #solve#28
    @ C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:899 [inlined]
 [13] solve(prob::NonlinearProblem{Vector{Float64}, true, SciMLBase.NullParameters, NonlinearFunction{true,SciMLBase.AutoSpecialize,…}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardNonlinearProblem}, args::NewtonRaphson{0, false, Val{:forward}, KrylovJL{typeof(gmres!), Int64, Tuple{}, Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, typeof(NonlinearSolve.DEFAULT_PRECS), true, nothing})
    @ DiffEqBase C:\Users\User\.julia\packages\DiffEqBase\JH4gt\src\solve.jl:889
 [14] top-level scope
    @ Untitled-1:12

The issue in the latter case seems to be from the jac_config field in the cache being empty https://github.com/SciML/NonlinearSolve.jl/blob/master/src/jacobian.jl#L27.

These errors also occur with TrustRegion. Seems to be fine with e.g. DynamicSS(TRBDF2(linsolve=KrylovJL_GMRES()))) from SteadyStateDiffEq.

[ChrisRackauckas]

I hadn't gotten to optimizing the Krylov setup yet since indeed they need a few things to be different from the direct methods. See #140 . Though I am a bit surprised to see them fail. But yes in general Krylov isn't ready here yet, but it's not too far off.

[ChrisRackauckas]

Handled in #203