Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Pre-allocate d and add factorization check #82

Merged
merged 2 commits into from
Nov 28, 2022
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 19 additions & 11 deletions src/CaNNOLeS.jl
Original file line number Diff line number Diff line change
Expand Up @@ -88,11 +88,13 @@ stats
"""
mutable struct CaNNOLeSSolver{T, V} <: AbstractOptimizationSolver
x::V
d::V
end

function CaNNOLeSSolver(nls::AbstractNLSModel{T, V}) where {T, V}
x = similar(nls.meta.x0)
return CaNNOLeSSolver{T, V}(x)
d = zeros(T, nls.meta.nvar + nls.nls_meta.nequ + nls.meta.ncon)
return CaNNOLeSSolver{T, V}(x, d)
end

function SolverCore.reset!(solver::CaNNOLeSSolver)
Expand Down Expand Up @@ -249,8 +251,9 @@ function SolverCore.solve!(
Jcx = ncon > 0 ? sparse(Jcx_rows, Jcx_cols, Jcx_vals, ncon, nvar) : spzeros(ncon, nvar)

r = copy(Fx)
dx = zeros(T, nvar)
dr = zeros(T, nequ)
d = solver.d
dx = view(d, 1:nvar)
dr = view(d, nvar .+ (1:nequ))
dλ = zeros(T, ncon)

Jxtr = Jx' * r
Expand Down Expand Up @@ -398,13 +401,15 @@ function SolverCore.solve!(
# on first time, μnew = μ⁺
rhs .= [dual; primal]
d, newton_success, ρ, ρold, nfacti =
newton_system(nvar, nequ, ncon, rhs, LDLT, ρold, params, linsolve)
newton_system!(d, nvar, nequ, ncon, rhs, LDLT, ρold, params, linsolve)
nfact += nfacti
nlinsolve += 1

if ρ > params[:ρmax] || any(isinf.(d)) || any(isnan.(d)) || fx ≥ T(1e60) # Error on hs70
if ρ > params[:ρmax] || !newton_success || any(isinf.(d)) || any(isnan.(d)) || fx ≥ T(1e60) # Error on hs70
internal_msg = if ρ > params[:ρmax]
"ρ → ∞"
elseif !newton_success
"Failure in Newton step computation"
elseif any(isinf.(d))
"d → ∞"
elseif any(isnan.(d))
Expand All @@ -413,11 +418,9 @@ function SolverCore.solve!(
"f → ∞"
end
broken = true
continue
break
end

dx .= d[1:nvar]
dr .= d[nvar .+ (1:nequ)]
dλ .= -d[nvar .+ nequ .+ (1:ncon)]
end # inner_iter != 1
### End of System solution
Expand Down Expand Up @@ -636,7 +639,7 @@ end
)

"""
newton_system(nvar, nequ, ncon, rhs, LDLT, ρold, params, linsolve)
newton_system!(d, nvar, nequ, ncon, rhs, LDLT, ρold, params, linsolve)

Compute an LDLt factorization of the (`nvar + nequ + ncon`)-square matrix for the Newton system contained in `LDLT`, i.e., `sparse(LDLT.rows, LDLT.cols, LDLT.vals, N, N)`, with the method `linsolve`.
If the factorization fails, a new factorization is attempted with an increased value for the regularization ρ as long as it is smaller than `params[:ρmax]`.
Expand All @@ -650,7 +653,8 @@ The factorization is then used to solve the linear system whose right-hand side
- `ρold`: the value of the regularization parameter used in the previous successful factorization, or 0 if this is the first one;
- `nfact`: the number of factorization attempts.
"""
function newton_system(
function newton_system!(
d::AbstractVector{T},
nvar::Integer,
nequ::Integer,
ncon::Integer,
Expand Down Expand Up @@ -722,7 +726,11 @@ function newton_system(
end
end

d, solve_success = solve(rhs, LDLT.factor)
if success
d, solve_success = solve!(rhs, LDLT.factor, d)
else
solve_success = false
end

return d, solve_success, ρ, ρold, nfact
end
Expand Down
10 changes: 5 additions & 5 deletions src/solver_types.jl
Original file line number Diff line number Diff line change
Expand Up @@ -12,8 +12,8 @@ if isdefined(HSL, :libhsl_ma57)
end

get_vals(LDLT::MA57Struct) = LDLT.factor.vals
function solve(rhs::AbstractVector, factor::Ma57)
d = ma57_solve(factor, -rhs)
function solve!(rhs::AbstractVector, factor::Ma57, d::AbstractVector)
d .= ma57_solve(factor, -rhs)
return d, true
end
else
Expand All @@ -35,8 +35,8 @@ end

get_vals(LDLT::LinearSolverStruct) = LDLT.vals

solve(::AbstractVector, factor::Nothing) = error("LDLt factorization failed.")
function solve(rhs::AbstractVector, factor::LDLFactorizations.LDLFactorization)
d = -(factor \ rhs)
solve!(::AbstractVector, factor::Nothing, ::AbstractVector) = error("LDLt factorization failed.")
function solve!(rhs::AbstractVector, factor::LDLFactorizations.LDLFactorization, d::AbstractVector)
d .= -(factor \ rhs)
return d, true
end