Sundials ERROR: UndefVarError: `LS` not defined

[briochemc]

Just tested solve(prob, KINSOL(;linear_solver=:GMRES)) on a simple problem and it threw:

ERROR: UndefVarError: `LS` not defined
Stacktrace:
  [1] ___kinsol(f::Function, y0::Vector{Float64}; userdata::Nothing, linear_solver::Symbol, jac_upper::Int64, jac_lower::Int64)
    @ Sundials ~/.julia/packages/Sundials/vGumE/src/simple.jl:15

[avik-pal]

Sundials.jl/src/simple.jl

Lines 70 to 77 in b01e0c3

	if linear_solver == :Dense
	A = Sundials.SUNDenseMatrix(length(y0), length(y0))
	LS = Sundials.SUNLinSol_Dense(y0, A)
	elseif linear_solver == :Band
	A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
	LS = Sundials.SUNLinSol_Band(y0, A)
	end
	flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A) true

seems like it doesn't have a case for GMRES. @ChrisRackauckas might be able to say more on this

[avik-pal]

On a sidenote, for most of our benchmarks we beat the nonlinear solver in Sundials by quite a margin.

[briochemc]

Sundials.jl/src/simple.jl

Lines 70 to 77 in b01e0c3

	if linear_solver == :Dense
	A = Sundials.SUNDenseMatrix(length(y0), length(y0))
	LS = Sundials.SUNLinSol_Dense(y0, A)
	elseif linear_solver == :Band
	A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
	LS = Sundials.SUNLinSol_Band(y0, A)
	end
	flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A) true

seems like it doesn't have a case for GMRES. @ChrisRackauckas might be able to say more on this

FWIW, I tried :KLU and it threw the same error.

[briochemc]

On a sidenote, for most of our benchmarks we beat the nonlinear solver in Sundials by quite a margin.

Do these benchmarks include KINSOL's GMRES or KLU?

[ChrisRackauckas]

They do not yet, but it's not clear why that would change anything. As a quasi-Newton algorithm it would see less optimizations in the GMRES case, so that's a worst case scenario for it. With KLU it's still a non-globalizing method. See:

https://docs.sciml.ai/SciMLBenchmarksOutput/dev/NonlinearProblem/nonlinear_solver_23_tests/

See, the point is that it doesn't have a globalizing part (line search, trust region, etc.) and makes heavy convergence assumptions in a quasi-Newton form. So it's issue isn't speed (though it's not the fastest) as much as convergence. It only successfully solves 15 of the 23 test cases, while some of the native NonlinearSolve.jl algorithms are both faster and converge on all 23. From all of the testing we've been able to do, we have not seen KINSOL as competitive at all, so it effectively only exists as a benchmarking tool and not something we'd ever recommend in practice.

[ChrisRackauckas]

But yes, if someone wants this, you just need to do

Sundials.jl/src/common_interface/solve.jl

Lines 263 to 295 in b01e0c3

	elseif LinearSolver == :Diagonal
	nojacobian = false
	flag = CVDiag(mem)
	_A = nothing
	_LS = nothing
	elseif LinearSolver == :GMRES
	LS = SUNLinSol_SPGMR(uvec, alg.prec_side, alg.krylov_dim)
	_A = nothing
	_LS = Sundials.LinSolHandle(LS, Sundials.SPGMR())
	elseif LinearSolver == :FGMRES
	LS = SUNLinSol_SPFGMR(uvec, alg.prec_side, alg.krylov_dim)
	_A = nothing
	_LS = LinSolHandle(LS, SPFGMR())
	elseif LinearSolver == :BCG
	LS = SUNLinSol_SPBCGS(uvec, alg.prec_side, alg.krylov_dim)
	_A = nothing
	_LS = LinSolHandle(LS, SPBCGS())
	elseif LinearSolver == :PCG
	LS = SUNLinSol_PCG(uvec, alg.prec_side, alg.krylov_dim)
	_A = nothing
	_LS = LinSolHandle(LS, PCG())
	elseif LinearSolver == :TFQMR
	LS = SUNLinSol_SPTFQMR(uvec, alg.prec_side, alg.krylov_dim)
	_A = nothing
	_LS = LinSolHandle(LS, PTFQMR())
	elseif LinearSolver == :KLU
	nojacobian = false
	nnz = length(SparseArrays.nonzeros(prob.f.jac_prototype))
	A = SUNSparseMatrix(length(uvec), length(uvec), nnz, CSC_MAT)
	LS = SUNLinSol_KLU(uvec, A)
	_A = MatrixHandle(A, SparseMatrix())
	_LS = LinSolHandle(LS, KLU())
	end

in the KINSOL setup code

Sundials.jl/src/simple.jl

Lines 70 to 77 in b01e0c3

	if linear_solver == :Dense
	A = Sundials.SUNDenseMatrix(length(y0), length(y0))
	LS = Sundials.SUNLinSol_Dense(y0, A)
	elseif linear_solver == :Band
	A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
	LS = Sundials.SUNLinSol_Band(y0, A)
	end
	flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A) true

. It would be easy to do, so it would be a nice contribution for someone to make at least so we can add GMRES and KLU to the benchmarks there, though since there isn't a globalizer in its method that wouldn't effect (improve) the convergence at all just performance in large-scale scenarios.

[jguterl]

Below is a very ugly draft implementation of gmres for kinsol with precondition. kinsol_prec is the entry point. using Sundials import Sundials: N_Vector, SUNMatrix, @checkflag, KIN_SUCCESS, UserFunctionAndData, PREC_RIGHT, PREC_NONE, KINCreate, KINSetUserData, KINInit, Handle, NVector using LinearAlgebra using SparseArrays using IncompleteLU abstract type AbstractKinsolData{J} end abstract type AbstractKinsolFunData{J,P} end abstract type AbstractKinsolPreconditioner{M} end abstract type AbstractPreconditionerSettings end struct ILU end struct AlgebraicMultiGrid end struct KinsolSettings linear_solver::Symbol y0::Vector jac_lower::Int64 jac_upper::Int64 prec_side::Int64 krylov_dim::Int64 print_level::Int64 end # struct KinsolFunData{F,D,J,P} <: AbstractKinsolFunData{J,P} # func::F # data::D # end struct KinsolData{J,P} jac::J preconditioner::P end # KinsolFunData(func::F, kinsol_data::KinsolData{J,P}) where {F,J,P} = KinsolFunData{F,D,J,P}(func, kinsol_data) # KinsolFunData(func::F, jac::J, precond::P) where {F,J<:Array,P} = KinsolFunData(func, KinsolData(jac, precond)) struct KinsolJacobian end ILUKinsolPreconditioner = AbstractKinsolPreconditioner{ILU} AlgebraicMultiGridKinsolPreconditioner = AbstractKinsolPreconditioner{AlgebraicMultiGrid} mutable struct KinsolPreconditioner{A,J,S,M} <: AbstractKinsolPreconditioner{M} P::A compute_jacobian::J settings::S method::M end @kwdef mutable struct AlgebraicMultiGridSettings <: AbstractPreconditionerSettings τ::Float64 = 1e-6 end @kwdef mutable struct ILUSettings <: AbstractPreconditionerSettings τ::Float64 = 1e-6 end function KinsolPreconditioner(jac, method::Symbol; kw...) if method == :ilu settings = ILUSettings(; kw...) P = compute_P_ilu(jac.jac, settings) method_ = ILU() KinsolPreconditioner(P, jac, settings, method_) elseif method == :multigrid settings = AlgebraicMultiGridSettings(; kw...) P = compute_P_AlgebraicMultiGrid(jac.jac, settings) method_ = AlgebraicMultiGrid() KinsolPreconditioner(P, jac, settings, method_) else error() end end function compute_P!(precond::T) where {T<:ILUKinsolPreconditioner} println("precond.settings.τ = $(precond.settings.τ)") precond.P = ilu(precond.compute_jacobian.jac, τ=precond.settings.τ) end compute_P_ilu(jac, settings) = ilu(jac, τ=settings.τ) compute_P_AlgebraicMultiGrid(jac, settings) = aspreconditioner( smoothed_aggregation(jac, presmoother=AlgebraicMultigrid.Jacobi(rand(size(jac, 1))), postsmoother=AlgebraicMultigrid.Jacobi(rand(size(jac, 1))))) psetup(u, uscale, fval, fscal, data::KinsolData{J,P}) where {J,P<: KinsolPreconditioner} = setup_P!(u, data.preconditioner) psolve(u, uscale, fval, fscal, v, data::KinsolData{J,P}) where {J,P<: KinsolPreconditioner} = solve!(data.preconditioner, u, uscale, fval, fscal, v) function _psetup(u::N_Vector, uscale::N_Vector, fval::N_Vector, fscale::N_Vector, userfun::UserFunctionAndData) psetup(convert(Vector, u), convert(Vector, uscale), convert(Vector, fval), convert(Vector, fscale), userfun.data) return KIN_SUCCESS end function _psolve(u::N_Vector, uscale::N_Vector, f::N_Vector, fscale::N_Vector, v::N_Vector, userfun::UserFunctionAndData ) psolve(convert(Vector, u), convert(Vector, uscale), convert(Vector, f), convert(Vector, fscale), convert(Vector, v), userfun.data) return KIN_SUCCESS end function compute_jac!(pc::KinsolPreconditioner{A,J,S,M}, u::Vector{Float64}) where {A,J,S,M} println("computing jacobian") pc.compute_jacobian(u) end import Sundials: kinsolfun # function kinsolfun(y::N_Vector, fy::N_Vector, userfun::KinsolFunData) # userfun[].func(convert(Vector, fy), convert(Vector, y), userfun[].data) # return KIN_SUCCESS # end function _kinsoljac( u::N_Vector, fu::N_Vector, J::SUNMatrix, userfun::UserFunctionAndData, tmp1::N_Vector, tmp2::N_Vector) userfun.data.preconditioner.compute_jacobian(convert(Vector, u)) copyto!(J, userfun.data.preconditioner.compute_jacobian.jac) return KIN_SUCCESS end function kinsol_prec(f, y0::Vector{Float64}; preconditioner::Union{Nothing,AbstractKinsolPreconditioner,Symbol}=nothing, jac::Union{KinsolJacobian,Nothing}=nothing, set_jacobian=false, linear_solver=:Band, ftol=nothing, jac_upper=0, jac_lower=0, su_scale=nothing, sf_scale=nothing, printlevel::Int64=1, krylov_dim=0, constraints=zeros(length(y0)), strategy=:none, max_iter=25 ) uvec = vec(deepcopy(y0)) userdata = KinsolData(jac, preconditioner) # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64}) # where `y` is the input vector, and `fy` is the result of the function # y0, Vector of initial values # return: the solution vector # --- initialize --- mem_ptr = KINCreate() (mem_ptr == C_NULL) && error("Failed to allocate KINSOL solver object") kmem = Handle(mem_ptr) # use the user_data field to pass a function # see: JuliaLang/julia#2554 userfun = UserFunctionAndData{typeof(f)}(f, userdata) function getcfun(::T) where {T} @cfunction(kinsolfun, Cint, (N_Vector, N_Vector, Ref{T})) end flag = @checkflag KINInit(kmem, getcfun(userfun), NVector(y0)) true flag = @checkflag KINSetUserData(kmem, userfun) true # this is where the F(U) is passed to kinsol through userfun pointer # --- preconditioner --- if preconditioner isa KinsolPreconditioner prec_side = PREC_RIGHT else prec_side = PREC_NONE end @Assert prec_side ∈ [PREC_NONE, PREC_RIGHT] # --- linear solver --- # J is template for Jacobian if linear_solver == :Dense J = Sundials.SUNDenseMatrix(neq, length(uvec)) LS = Sundials.SUNLinSol_Dense(uvec, L) elseif linear_solver == :LapackDense J = Sundials.SUNDenseMatrix(length(uvec), length(uvec)) LS = Sundials.SUNLinSol_LapackDense(uvec, J) elseif linear_solver == :Band J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower) LS = Sundials.SUNLinSol_Band(uvec, J) elseif linear_solver == :LapackBand J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower) LS = Sundials.SUNLinSol_LapackBand(uvec, J) elseif linear_solver == :GMRES LS = Sundials.SUNLinSol_SPGMR(uvec, prec_side, krylov_dim) J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower) elseif linear_solver == :FGMRES LS = Sundials.SUNLinSol_SPFGMR(uvec, prec_side, krylov_dim) J = Sundials.SUNBandMatrix(length(uvec), jac_upper, jac_lower) elseif linear_solver == :KLU nnz = length(SparseArrays.nonzeros(preconditioner.compute_jacobian.jac)) J = Sundials.SUNSparseMatrix(length(uvec), length(uvec), nnz, Sundials. CSC_MAT) LS = Sundials.SUNLinSol_KLU(uvec, J) else error("linear solver available for Kinsol are: [:Dense,:LapackDense,:Band, :GMRES, :FGMRES, :KLU]") end flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, J === nothing ? C_NULL : J) true J === nothing ? jacobian = false : jacobian = true # --- jacobian --- if preconditioner isa KinsolPreconditioner && set_jacobian println("Kinsol: setting up Jacobian ...") function getjacfun(::T) where {T} @cfunction(_kinsoljac, Cint, ( N_Vector, N_Vector, SUNMatrix, Ref{T}, N_Vector, N_Vector)) end jac = getjacfun(userfun) flag = @checkflag Sundials.KINSetJacFn(kmem, jac) true end function getpsolvefun(::T) where {T} @cfunction(_psolve, Cint, ( N_Vector, N_Vector, N_Vector, N_Vector, N_Vector, Ref{T})) end function getpsetupfun(::T) where {T} @cfunction(_psetup, Cint, ( N_Vector, N_Vector, N_Vector, N_Vector, Ref{T})) end psolvefun = getpsolvefun(userfun) psetupfun = getpsetupfun(userfun) # --- preconditionner setting --- if preconditioner isa AbstractKinsolPreconditioner && linear_solver != :KLU println("KINSOL: setting preconditioner...") Sundials.KINSetPreconditioner(kmem, psetupfun, psolvefun) end ## Solve problem if isnothing(su_scale) su_scale = ones(length(y0)) end if isnothing(sf_scale) sf_scale = ones(length(y0)) end @Assert length(sf_scale) == length(y0) @Assert length(su_scale) == length(y0) if strategy == :none strategy = Sundials.KIN_NONE elseif strategy == :linesearch strategy = Sundials.KIN_LINESEARCH elseif strategy == :picard strategy = Sundials.KIN_PICARD else error("strategy $strategy not available for kinsol. Available strategies:$([:none, :lineasearch, :picard])") end # jacobian Sundials.KINSetPrintLevel(kmem, printlevel) if ftol isa Number Sundials.KINSetFuncNormTol(kmem, ftol) end Sundials.KINSetNumMaxIters(kmem, max_iter) Sundials.KINSetConstraints(kmem, constraints) etachoice = Sundials.KIN_ETACHOICE2 return kmem end function run_kinsol(kmem, y, strategy::Union{Symbol,Int32}, su_scale:: Union{Nothing,Vector{Float64}}, sf_scale::Union{Nothing,Vector{Float64}}) if isnothing(su_scale) su_scale_ = ones(length(y)) else su_scale_ = su_scale end if isnothing(sf_scale) sf_scale = ones(length(y)) end y_ = deepcopy(y) #Sundials.KINSetEtaForm(kmem, etachoice) #println("su_scale:", su_scale) #println("sf_scale:", sf_scale) flag = @checkflag Sundials.KINSol(kmem, y_, kinsol_strategy(strategy), su_scale_, sf_scale) true # Sundials.SUNLinSolFree_Dense(LS) # Sundials.SUNMatDestroy_Dense(J) return (y_, flag) end kinsol_strategy(strategy::Int32) = strategy function kinsol_strategy(strategy::Symbol) if strategy == :none strategy = Sundials.KIN_NONE elseif strategy == :linesearch strategy = Sundials.KIN_LINESEARCH elseif strategy == :picard strategy = Sundials.KIN_PICARD else error("strategy $strategy not available for kinsol. Available strategies:$([:none, :lineasearch, :picard])") end return strategy end struct KinSOLPreconditioner{MP,J,F,S,M} <: AbstractKinsolPreconditioner{M} P::MP jac::J evaluate_jac!::F settings::S method::M end function kinsolfun(y::N_Vector, fy::N_Vector, userfun::UserFunctionAndData) userfun.func(convert(Vector, fy), convert(Vector, y), userfun.data) return Sundials.KIN_SUCCESS end function setup_P!(u, precond) #println("setup_P") # Build preconditioner compute_jac!(precond, u) compute_P!(precond) end function solve!(precond::T, u, uscale, fval, fscal, v) where {T<: ILUKinsolPreconditioner} #println("solve Px=c") #p = fscal .* v ldiv!(precond.P, v) #v .= p ./ fscal end function solve!(precond::T, u, uscale, fval, fscal, v) where {T<: AlgebraicMultiGridKinsolPreconditioner} #println("solve Px=c") #p = fscal .* v ldiv!(precond.P, v) #v .= p ./ fscal end using AlgebraicMultigrid function compute_P!(precond::T) where {T<: AlgebraicMultiGridKinsolPreconditioner} precond.P = aspreconditioner(smoothed_aggregation(precond.compute_jacobian.jac, presmoother=AlgebraicMultigrid.Jacobi(rand(size(precond.compute_jacobian.jac, 1))), postsmoother=AlgebraicMultigrid.Jacobi(rand(size(precond. compute_jacobian.jac, 1))))) end # function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata) # if newW === nothing || newW # Pl = aspreconditioner(ruge_stuben(convert(AbstractMatrix, W))) # else # Pl = Plprev # end # Pl, nothing # end ############# ENTRY POINT (out_ks, flag) = kinsol_prec((f, u, data) -> f!(f, u), U; linear_solver=:GMRES, printlevel=0, krylov_dim=floor(Int64, sqrt(length(U))), strategy=:none, max_iter=50)

…

On Wed, Nov 1, 2023 at 10:33 PM Benoit Pasquier ***@***.***> wrote: Just tested solve(prob, KINSOL(;linear_solver=:GMRES)) on a simple problem and it threw: ERROR: UndefVarError: `LS` not defined Stacktrace: [1] ___kinsol(f::Function, y0::Vector{Float64}; userdata::Nothing, linear_solver::Symbol, jac_upper::Int64, jac_lower::Int64) @ Sundials ~/.julia/packages/Sundials/vGumE/src/simple.jl:15 — Reply to this email directly, view it on GitHub <#431>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AEESMZDWFJHL4VTU6UZRW2TYCMV7TAVCNFSM6AAAAAA62KHFJWVHI2DSMVQWIX3LMV43ASLTON2WKOZRHE3TGNJTHEYTKNI> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>