JuliaNLSolvers / NLSolvers.jl

No bells and whistles foundation of Optim.jl
https://julianlsolvers.github.io/NLSolvers.jl/dev/
MIT License
27 stars 8 forks source link

Setup with LinearSolve.jl #24

Open ChrisRackauckas opened 2 years ago

ChrisRackauckas commented 2 years ago

I'm trying to get more libraries using the LinearSolve.jl interface

http://linearsolve.sciml.ai/dev/

It's pretty complete now. The suggested form is that you use the caching interface so that it can store factorizations and iterative solver inits (http://linearsolve.sciml.ai/dev/tutorials/caching_interface/). The suggested algorithm form is to just have the user pass the type, like Newton(; linsolve=KLU()) which would then switch over to using SuiteSparse.KLU and accelerate sparse matrix factorizations for unstructured cases over UMFPACK. http://linearsolve.sciml.ai/dev/solvers/solvers/#SuiteSparse.jl

longemen3000 commented 2 years ago

one way to do this is to delay the constructor of the linsolve struct until inside the solve function, something like:

struct LinearSolveCache{T}
  linsolve::T
end

function (cache::LinearSolveCache)(d, B, g)
  LinearSolve.set_a(cache.linsolve.cache,B)
  LinearSolve.set_b(cache.linsolve.cache,g)
  solve(cache.linsolve)
  d .= cache.linsolve.u
end

#does it work with static arrays?
function (cache::LinearSolveCache)(B, g)
  prob = LinearProblem(B,g)
  solve(prob,cache.linsolve)
  prob.u
end

init_linsolve(linsolve::T,A,b) where T <: Function = linsolve
init_linsolve(linsolve:::LinearSolve.SciMLLinearSolveAlgorithm) = LinearSolveCache(init(LinearProblem(A,b),linsolve))

and then add those inits in the newton solver for NeqProblem and OptimizationProblem (the last one seems harder, but it could be added to the objvars result? what do you think @pkofod ?