diff --git a/scripts/common.jl b/scripts/common.jl index 916c0f35..71b6d527 100644 --- a/scripts/common.jl +++ b/scripts/common.jl @@ -5,6 +5,7 @@ using CodecBzip2 using HSL using NLPModels using SparseArrays +using Quadmath import QuadraticModels: SparseMatrixCOO @@ -43,6 +44,49 @@ function _scale_coo!(A, Dr, Dc) end end +function presolve(qp; max_iter=3) + pqp = qp + for i in 1:max_iter + println(NLPModels.get_nvar(pqp)) + pqp = MadIPM.presolve_qp(pqp)[1] + end + return pqp +end + +function _scale_qp(qp::QuadraticModel, Dr, Dc) + Hs = copy(qp.data.H) + As = copy(qp.data.A) + _scale_coo!(Hs, Dc, Dc) + _scale_coo!(As, Dr, Dc) + + data = QuadraticModels.QPData( + qp.data.c0, + qp.data.c ./ Dc, + Hs, + As, + ) + + return QuadraticModel( + NLPModelMeta( + qp.meta.nvar; + ncon=qp.meta.ncon, + lvar=qp.meta.lvar .* Dc, + uvar=qp.meta.uvar .* Dc, + lcon=qp.meta.lcon ./ Dr, + ucon=qp.meta.ucon ./ Dr, + x0=qp.meta.x0 .* Dc, + y0=qp.meta.y0 ./ Dr, + nnzj=qp.meta.nnzj, + lin_nnzj=qp.meta.nnzj, + lin=qp.meta.lin, + nnzh=qp.meta.nnzh, + minimize=qp.meta.minimize, + ), + Counters(), + data, + ) +end + """ scale_qp(qp::QuadraticModel) @@ -65,29 +109,26 @@ function scale_qp(qp::QuadraticModel) A_csc = sparse(A.rows, A.cols, A.vals, m, n) Dr, Dc = HSL.mc77(A_csc, 0) - Hs = copy(qp.data.H) - As = copy(qp.data.A) - _scale_coo!(Hs, Dc, Dc) - _scale_coo!(As, Dr, Dc) + return _scale_qp(qp, Dr, Dc) +end +function convert_qp(qp::QuadraticModel, T::Type{<:AbstractFloat}) data = QuadraticModels.QPData( - qp.data.c0, - qp.data.c ./ Dc, - qp.data.v, - Hs, - As, + T(qp.data.c0), + convert(Vector{T}, qp.data.c), + convert(SparseMatrixCOO{T, Int}, qp.data.H), + convert(SparseMatrixCOO{T, Int}, qp.data.A), ) - return QuadraticModel( NLPModelMeta( qp.meta.nvar; ncon=qp.meta.ncon, - lvar=qp.meta.lvar .* Dc, - uvar=qp.meta.uvar .* Dc, - lcon=qp.meta.lcon ./ Dr, - ucon=qp.meta.ucon ./ Dr, - x0=qp.meta.x0 .* Dc, - y0=qp.meta.y0 ./ Dr, + lvar=convert(Vector{T}, qp.meta.lvar), + uvar=convert(Vector{T}, qp.meta.uvar), + lcon=convert(Vector{T}, qp.meta.lcon), + ucon=convert(Vector{T}, qp.meta.ucon), + x0=convert(Vector{T}, qp.meta.x0), + y0=convert(Vector{T}, qp.meta.y0), nnzj=qp.meta.nnzj, lin_nnzj=qp.meta.nnzj, lin=qp.meta.lin, @@ -99,3 +140,4 @@ function scale_qp(qp::QuadraticModel) ) end + diff --git a/src/kernels.jl b/src/kernels.jl index 88cd5687..74f3f66e 100644 --- a/src/kernels.jl +++ b/src/kernels.jl @@ -1,6 +1,6 @@ -function set_initial_primal_rhs!(solver::MadNLP.AbstractMadNLPSolver) +function set_initial_primal_rhs!(solver::MadNLP.AbstractMadNLPSolver{T}) where T p = solver.p - fill!(full(p), 0.0) + fill!(full(p), zero(T)) py = MadNLP.dual(p) b = solver.c @@ -8,9 +8,9 @@ function set_initial_primal_rhs!(solver::MadNLP.AbstractMadNLPSolver) return end -function set_initial_dual_rhs!(solver::MadNLP.AbstractMadNLPSolver) +function set_initial_dual_rhs!(solver::MadNLP.AbstractMadNLPSolver{T}) where T p = solver.p - fill!(full(p), 0.0) + fill!(full(p), zero(T)) px = MadNLP.primal(p) c = MadNLP.primal(solver.f) @@ -18,7 +18,7 @@ function set_initial_dual_rhs!(solver::MadNLP.AbstractMadNLPSolver) return end -function set_predictive_rhs!(solver::MadNLP.AbstractMadNLPSolver, kkt::MadNLP.AbstractKKTSystem) +function set_predictive_rhs!(solver::MadNLP.AbstractMadNLPSolver{T}, kkt::MadNLP.AbstractKKTSystem{T}) where T # RHS px = MadNLP.primal(solver.p) py = MadNLP.dual(solver.p) @@ -31,7 +31,7 @@ function set_predictive_rhs!(solver::MadNLP.AbstractMadNLPSolver, kkt::MadNLP.Ab # Constraint c = solver.c - fill!(MadNLP.full(solver.p), 0.0) + fill!(MadNLP.full(solver.p), zero(T)) px .= .-f .+ zl .- zu .- solver.jacl py .= .-c @@ -40,7 +40,7 @@ function set_predictive_rhs!(solver::MadNLP.AbstractMadNLPSolver, kkt::MadNLP.Ab return end -function set_correction_rhs!(solver::MadNLP.AbstractMadNLPSolver, kkt::MadNLP.AbstractKKTSystem, mu::Float64, correction_lb::AbstractVector{Float64}, correction_ub::AbstractVector{Float64}, ind_lb, ind_ub) +function set_correction_rhs!(solver::MadNLP.AbstractMadNLPSolver, kkt::MadNLP.AbstractKKTSystem, mu::T, correction_lb::AbstractVector{T}, correction_ub::AbstractVector{T}, ind_lb, ind_ub) where T px = MadNLP.primal(solver.p) py = MadNLP.dual(solver.p) pzl = MadNLP.dual_lb(solver.p) @@ -72,10 +72,10 @@ end # Gondzio's multi-correction scheme function set_extra_correction!( - solver::MadNLP.AbstractMadNLPSolver, + solver::MadNLP.AbstractMadNLPSolver{T}, correction_lb, correction_ub, alpha_p, alpha_d, βmin, βmax, μ, -) +) where T dlb = MadNLP.dual_lb(solver.d) dub = MadNLP.dual_ub(solver.d) tmin, tmax = βmin * μ, βmax * μ @@ -91,7 +91,7 @@ function set_extra_correction!( elseif v > tmax tmax - v else - 0.0 + zero(T) end corr - δ end, @@ -110,7 +110,7 @@ function set_extra_correction!( elseif v > tmax tmax - v else - 0.0 + zero(T) end corr + δ end, @@ -174,10 +174,10 @@ end Barrier =# -function get_complementarity_measure(solver::MadNLP.AbstractMadNLPSolver) +function get_complementarity_measure(solver::MadNLP.AbstractMadNLPSolver{T}) where T m1, m2 = length(solver.x_lr), length(solver.x_ur) if m1 + m2 == 0 - return 0.0 + return zero(T) else inf_compl_l = mapreduce( (x_lr, xl_r, zl_r) -> (x_lr - xl_r) * zl_r, @@ -195,10 +195,10 @@ function get_complementarity_measure(solver::MadNLP.AbstractMadNLPSolver) end end -function get_affine_complementarity_measure(solver::MadNLP.AbstractMadNLPSolver, alpha_p, alpha_d) +function get_affine_complementarity_measure(solver::MadNLP.AbstractMadNLPSolver{T}, alpha_p, alpha_d) where T m1, m2 = length(solver.x_lr), length(solver.x_ur) if m1 + m2 == 0 - return 0.0 + return zero(T) else dzlb = MadNLP.dual_lb(solver.d) dzub = MadNLP.dual_ub(solver.d) @@ -327,9 +327,9 @@ end # Implement Mehrotra's heuristic to compute the step : see Procedure GTSF (Exhibit 6.1) in # "On The Implementation Of A Primal-Dual Interior Point Method" -function update_step!(rule::MehrotraAdaptiveStep, solver) - gamma_a = 1.0 / (1.0 - rule.gamma_f) - tau = 1.0 +function update_step!(rule::MehrotraAdaptiveStep, solver::MPCSolver{T}) where T + gamma_a = one(T) / (one(T) - rule.gamma_f) + tau = one(T) d_zl = MadNLP.dual_lb(solver.d) d_zu = MadNLP.dual_ub(solver.d) @@ -349,11 +349,11 @@ function update_step!(rule::MehrotraAdaptiveStep, solver) mu_full = get_affine_complementarity_measure(solver, max_alpha_p, max_alpha_d) mu_full /= gamma_a - alpha_p, alpha_d = 1.0, 1.0 + alpha_p, alpha_d = one(T), one(T) # Require CUDA.jl as a dependency of MadIPM.jl # CUDA.@allowscalar begin - if max_alpha_p < 1.0 + if max_alpha_p < one(T) if alpha_xl <= alpha_xu tmp = mu_full / (solver.zl_r[i_xl] + max_alpha_d * d_zl[i_xl]) alpha_p = (solver.x_lr[i_xl] - solver.xl_r[i_xl] - tmp) / (-solver.dx_lr[i_xl]) @@ -362,7 +362,7 @@ function update_step!(rule::MehrotraAdaptiveStep, solver) alpha_p = (solver.xu_r[i_xu] - solver.x_ur[i_xu] - tmp) / (solver.dx_ur[i_xu]) end end - if max_alpha_d < 1.0 + if max_alpha_d < one(T) if alpha_zl <= alpha_zu tmp = mu_full / (solver.x_lr[i_zl] + max_alpha_p * solver.dx_lr[i_zl] - solver.xl_r[i_zl]) alpha_d = -(solver.zl_r[i_zl] - tmp) / d_zl[i_zl] @@ -382,20 +382,20 @@ end Regularization =# -function init_regularization!(solver::MPCSolver, ::NoRegularization) - solver.del_w = 1.0 - solver.del_c = 0.0 +function init_regularization!(solver::MPCSolver{T}, ::NoRegularization) where T + solver.del_w = one(T) + solver.del_c = zero(T) return end -function update_regularization!(solver::MPCSolver, ::NoRegularization) - solver.del_w = 0.0 - solver.del_c = 0.0 +function update_regularization!(solver::MPCSolver{T}, ::NoRegularization) where T + solver.del_w = zero(T) + solver.del_c = zero(T) return end -function init_regularization!(solver::MPCSolver, reg::FixedRegularization) - solver.del_w = 1.0 +function init_regularization!(solver::MPCSolver{T}, reg::FixedRegularization) where T + solver.del_w = one(T) solver.del_c = reg.delta_d return end @@ -406,16 +406,16 @@ function update_regularization!(solver::MPCSolver, reg::FixedRegularization) return end -function init_regularization!(solver::MPCSolver, reg::AdaptiveRegularization) - solver.del_w = 1.0 +function init_regularization!(solver::MPCSolver{T}, reg::AdaptiveRegularization) where T + solver.del_w = one(T) solver.del_c = reg.delta_d return end -function update_regularization!(solver::MPCSolver, reg::AdaptiveRegularization) - reg.delta_p = max(reg.delta_p / 10.0, reg.delta_min) +function update_regularization!(solver::MPCSolver{T}, reg::AdaptiveRegularization) where T + reg.delta_p = max(reg.delta_p / T(10.0), reg.delta_min) # Dual regularization is negative! - reg.delta_d = min(reg.delta_d / 10.0, -reg.delta_min) + reg.delta_d = min(reg.delta_d / T(10.0), -reg.delta_min) solver.del_w = reg.delta_p solver.del_c = reg.delta_d return @@ -437,7 +437,7 @@ function dual_objective(solver::MPCSolver) return dobj end -function get_optimality_gap(solver::MPCSolver) +function get_optimality_gap(solver::MPCSolver{T}) where T return MadNLP.get_inf_compl( solver.x_lr, solver.xl_r, @@ -445,7 +445,7 @@ function get_optimality_gap(solver::MPCSolver) solver.xu_r, solver.x_ur, solver.zu_r, - 0., - 1.0, + zero(T), + one(T), ) end diff --git a/src/linear_solver.jl b/src/linear_solver.jl index a6bc5e8f..454f391d 100644 --- a/src/linear_solver.jl +++ b/src/linear_solver.jl @@ -3,7 +3,7 @@ Interface to direct solver for solving KKT system =# -function factorize_regularized_system!(solver) +function factorize_regularized_system!(solver::MPCSolver{T}) where T max_trials = 3 for ntrial in 1:max_trials set_aug_diagonal_reg!(solver.kkt, solver) @@ -11,8 +11,8 @@ function factorize_regularized_system!(solver) if is_factorized(solver.kkt.linear_solver) break end - solver.del_w *= 100.0 - solver.del_c *= 100.0 + solver.del_w *= T(100) + solver.del_c *= T(100) end end @@ -27,6 +27,7 @@ function solve_system!( # Check residual w = solver._w1 + copyto!(MadNLP.full(w), MadNLP.full(p)) mul!(w, solver.kkt, d, -one(T), one(T)) norm_w = norm(MadNLP.full(w), Inf) diff --git a/src/solver.jl b/src/solver.jl index 8d8d32cc..718439e7 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -242,14 +242,14 @@ function mehrotra_correction_direction!(solver) return end -function gondzio_correction_direction!(solver) +function gondzio_correction_direction!(solver::MPCSolver{T}) where T solver.opt.max_ncorr ≤ 0 && return - δ = 0.1 - γ = 0.1 - βmin = 0.1 - βmax = 10.0 - tau = 0.995 + δ = T(0.1) + γ = T(0.1) + βmin = T(0.1) + βmax = T(10.0) + tau = T(0.995) # Load buffer for descent direction. Δp = solver._w2.values @@ -258,8 +258,8 @@ function gondzio_correction_direction!(solver) for ncorr in 1:solver.opt.max_ncorr # Enlarge step sizes in primal and dual spaces. - tilde_alpha_p = min(alpha_p + δ, 1.0) - tilde_alpha_d = min(alpha_d + δ, 1.0) + tilde_alpha_p = min(alpha_p + δ, one(T)) + tilde_alpha_d = min(alpha_d + δ, one(T)) # Apply Mehrotra's heuristic for centering parameter mu. ga = get_affine_complementarity_measure(solver, tilde_alpha_p, tilde_alpha_d) g = solver.mu_curr @@ -285,7 +285,7 @@ function gondzio_correction_direction!(solver) hat_alpha_p, hat_alpha_d = get_fraction_to_boundary_step(solver, tau) # Stop extra correction if the stepsize does not increase sufficiently - if (hat_alpha_p < 1.005 * alpha_p) || (hat_alpha_d < 1.005 * alpha_d) + if (hat_alpha_p < T(1.005) * alpha_p) || (hat_alpha_d < T(1.005) * alpha_d) copyto!(solver.d.values, Δp) break else diff --git a/src/utils.jl b/src/utils.jl index d8773830..7c68a557 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -66,8 +66,8 @@ end Options =# -@kwdef mutable struct IPMOptions <: MadNLP.AbstractOptions - tol::Float64 +@kwdef mutable struct IPMOptions{T} <: MadNLP.AbstractOptions + tol::T kkt_system::Type linear_solver::Type # Output options @@ -77,18 +77,18 @@ end rethrow_error::Bool = false # Termination options max_iter::Int = 3000 - max_wall_time::Float64 = 1e6 - divergence_tol::Float64 = 1e4 - divergence_scale::Float64 = 10.0 - kappa_d::Float64 = 1e-5 + max_wall_time::T = 1e6 + divergence_tol::T = 1e4 + divergence_scale::T = 10.0 + kappa_d::T = 1e-5 fixed_variable_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxBound : MadNLP.MakeParameter equality_treatment::Type = kkt_system <: MadNLP.SparseCondensedKKTSystem ? MadNLP.RelaxEquality : MadNLP.EnforceEquality # initialization options scaling::Bool = true - nlp_scaling_max_gradient::Float64 = 100.0 - bound_push::Float64 = 1e-2 - bound_fac::Float64 = 1e-2 - bound_relax_factor::Float64 = 1e-12 + nlp_scaling_max_gradient::T = 100.0 + bound_push::T = 1e-2 + bound_fac::T = 1e-2 + bound_relax_factor::T = 1e-12 # Regularization regularization::AbstractRegularization = FixedRegularization(1e-10, 1e-10) # Step @@ -96,13 +96,13 @@ end # Barrier barrier_update::AbstractBarrierUpdate = Mehrotra() max_ncorr::Int = 0 - s_max::Float64 = 100.0 - mu_init::Float64 = 1e-1 - mu_min::Float64 = 1e-12 - mu_superlinear_decrease_power::Float64 = 1.5 - tau_min::Float64 = 0.99 + s_max::T = 100.0 + mu_init::T = 1e-1 + mu_min::T = 1e-12 + mu_superlinear_decrease_power::T = 1.5 + tau_min::T = 0.99 # Linear solve - tol_linear_solve::Float64 = 1e-8 + tol_linear_solve::T = 1e-8 check_residual::Bool = false end @@ -111,9 +111,9 @@ function IPMOptions( nlp::NLPModels.AbstractNLPModel{T}; kkt_system = MadNLP.SparseKKTSystem, linear_solver = MadNLP.default_sparse_solver(nlp), - tol = 1e-8, + tol = T(1e-8), ) where T - return IPMOptions( + return IPMOptions{T}( tol = tol, kkt_system = kkt_system, linear_solver = linear_solver, @@ -467,6 +467,7 @@ function standard_form_qp(qp::QuadraticModels.QuadraticModel) ucon_[m + k] = xu[k] end + c_ = [qp.data.c; zeros(ns + nw)] lvar_ = [lvar; lcon[ind_ineq]; zeros(nw)] uvar_ = [uvar; ucon[ind_ineq]; fill(Inf, nw)] # The upper bounds in range constraints have been moved in a separate constraint