From 497eecf140523f78bf696c7b74e469efd6fc391b Mon Sep 17 00:00:00 2001 From: Neven Sajko Date: Sun, 9 Oct 2022 00:33:57 +0200 Subject: [PATCH] [perf] Fuse more vectorized operations (#127) --- src/IPM/HSD/HSD.jl | 12 +++++------- src/IPM/HSD/step.jl | 26 ++++++++++++-------------- src/IPM/MPC/MPC.jl | 24 +++++++++++------------- src/IPM/MPC/step.jl | 14 ++++++-------- src/IPM/ipmdata.jl | 4 ++-- 5 files changed, 36 insertions(+), 44 deletions(-) diff --git a/src/IPM/HSD/HSD.jl b/src/IPM/HSD/HSD.jl index 37efe7a4..e2681a11 100644 --- a/src/IPM/HSD/HSD.jl +++ b/src/IPM/HSD/HSD.jl @@ -88,14 +88,12 @@ function compute_residuals!(hsd::HSD{T} # Lower-bound residual # rl_j = τ*l_j - (x_j - xl_j) if l_j ∈ R # = 0 if l_j = -∞ - @. res.rl = - pt.x + pt.xl + pt.τ * dat.l - res.rl .*= dat.lflag + @. res.rl = (- pt.x + pt.xl + pt.τ * dat.l) * dat.lflag # Upper-bound residual # ru_j = τ*u_j - (x_j + xu_j) if u_j ∈ R # = 0 if u_j = +∞ - @. res.ru = - pt.x - pt.xu + pt.τ * dat.u - res.ru .*= dat.uflag + @. res.ru = (- pt.x - pt.xu + pt.τ * dat.u) * dat.uflag # Dual residual # rd = t*c - A'y - zl + zu @@ -173,8 +171,8 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher # Check for infeasibility certificates if max( norm(dat.A * pt.x, Inf), - norm((pt.x - pt.xl) .* dat.lflag, Inf), - norm((pt.x + pt.xu) .* dat.uflag, Inf) + norm((pt.x .- pt.xl) .* dat.lflag, Inf), + norm((pt.x .+ pt.xu) .* dat.uflag, Inf) ) * (norm(dat.c, Inf) / max(1, norm(dat.b, Inf))) < - ϵi * dot(dat.c, pt.x) # Dual infeasible, i.e., primal unbounded hsd.primal_status = Sln_InfeasibilityCertificate @@ -182,7 +180,7 @@ function update_solver_status!(hsd::HSD{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher return nothing end - δ = dat.A' * pt.y + (pt.zl .* dat.lflag) - (pt.zu .* dat.uflag) + δ = dat.A' * pt.y .+ (pt.zl .* dat.lflag) .- (pt.zu .* dat.uflag) if norm(δ, Inf) * max( norm(dat.l .* dat.lflag, Inf), norm(dat.u .* dat.uflag, Inf), diff --git a/src/IPM/HSD/step.jl b/src/IPM/HSD/step.jl index 5b48a7bf..cc7b7abc 100644 --- a/src/IPM/HSD/step.jl +++ b/src/IPM/HSD/step.jl @@ -58,7 +58,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst # Compute hx, hy, hz from first augmented system solve hx = tzeros(Tv, n) hy = tzeros(Tv, m) - ξ_ = dat.c - ((pt.zl ./ pt.xl) .* dat.l) .* dat.lflag - ((pt.zu ./ pt.xu) .* dat.u) .* dat.uflag + ξ_ = @. (dat.c - ((pt.zl / pt.xl) * dat.l) * dat.lflag - ((pt.zu / pt.xu) * dat.u) * dat.uflag) KKT.solve!(hx, hy, hsd.kkt, dat.b, ξ_) # Recover h0 = ρg + κ / τ - c'hx + b'hy - u'hz @@ -67,7 +67,7 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst h0 = ( dot(dat.l .* dat.lflag, (dat.l .* θl) .* dat.lflag) + dot(dat.u .* dat.uflag, (dat.u .* θu) .* dat.uflag) - - dot(c + (θl .* dat.l) .* dat.lflag + (θu .* dat.u) .* dat.uflag, hx) + - dot((@. (c + (θl * dat.l) * dat.lflag + (θu * dat.u) * dat.uflag)), hx) + dot(b, hy) + pt.κ / pt.τ + hsd.regG @@ -77,9 +77,9 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst @timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0, # Right-hand side of Newton system res.rp, res.rl, res.ru, res.rd, res.rg, - -(pt.xl .* pt.zl) .* dat.lflag, - -(pt.xu .* pt.zu) .* dat.uflag, - -pt.τ * pt.κ + .-(pt.xl .* pt.zl) .* dat.lflag, + .-(pt.xu .* pt.zu) .* dat.uflag, + .-pt.τ * pt.κ ) # Step length for affine-scaling direction @@ -91,8 +91,8 @@ function compute_step!(hsd::HSD{T, Tv}, params::IPMOptions{T}) where{T, Tv<:Abst @timeit hsd.timer "Newton" solve_newton_system!(Δ, hsd, hx, hy, h0, # Right-hand side of Newton system η .* res.rp, η .* res.rl, η .* res.ru, η .* res.rd, η * res.rg, - (-pt.xl .* pt.zl .+ γ * pt.μ .- Δ.xl .* Δ.zl) .* dat.lflag, - (-pt.xu .* pt.zu .+ γ * pt.μ .- Δ.xu .* Δ.zu) .* dat.uflag, + (.-pt.xl .* pt.zl .+ γ * pt.μ .- Δ.xl .* Δ.zl) .* dat.lflag, + (.-pt.xu .* pt.zu .+ γ * pt.μ .- Δ.xu .* Δ.zu) .* dat.uflag, -pt.τ * pt.κ + γ * pt.μ - Δ.τ * Δ.κ ) α = max_step_length(pt, Δ) @@ -222,9 +222,9 @@ function solve_newton_system!(Δ::Point{T, Tv}, @timeit hsd.timer "Δτ" Δ.τ = ( ξg_ - + dot(dat.c - + ((pt.zl ./ pt.xl) .* dat.l) .* dat.lflag - + ((pt.zu ./ pt.xu) .* dat.u) .* dat.uflag + + dot((@. (dat.c + + ((pt.zl / pt.xl) * dat.l) * dat.lflag + + ((pt.zu / pt.xu) * dat.u) * dat.uflag)) , Δ.x) - dot(dat.b, Δ.y) ) / h0 @@ -236,12 +236,10 @@ function solve_newton_system!(Δ::Point{T, Tv}, # III. Recover Δxl, Δxu @timeit hsd.timer "Δxl" begin - @. Δ.xl = -ξl + Δ.x - Δ.τ .* (dat.l .* dat.lflag) - Δ.xl .*= dat.lflag + @. Δ.xl = (-ξl + Δ.x - Δ.τ .* (dat.l .* dat.lflag)) * dat.lflag end @timeit hsd.timer "Δxu" begin - @. Δ.xu = ξu - Δ.x + Δ.τ .* (dat.u .* dat.uflag) - Δ.xu .*= dat.uflag + @. Δ.xu = ( ξu - Δ.x + Δ.τ .* (dat.u .* dat.uflag)) * dat.uflag end # IV. Recover Δzl, Δzu diff --git a/src/IPM/MPC/MPC.jl b/src/IPM/MPC/MPC.jl index b2ccefd5..55e46c22 100644 --- a/src/IPM/MPC/MPC.jl +++ b/src/IPM/MPC/MPC.jl @@ -111,14 +111,12 @@ function compute_residuals!(mpc::MPC{T}) where{T} # Lower-bound residual # rl_j = l_j - (x_j - xl_j) if l_j ∈ R # = 0 if l_j = -∞ - @. res.rl = (dat.l + pt.xl) - pt.x - res.rl .*= dat.lflag + @. res.rl = ((dat.l + pt.xl) - pt.x) * dat.lflag # Upper-bound residual # ru_j = u_j - (x_j + xu_j) if u_j ∈ R # = 0 if u_j = +∞ - @. res.ru = dat.u - (pt.x + pt.xu) - res.ru .*= dat.uflag + @. res.ru = (dat.u - (pt.x + pt.xu)) * dat.uflag # Dual residual # rd = c - (A'y + zl - zu) @@ -188,8 +186,8 @@ function update_solver_status!(mpc::MPC{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher # Check for infeasibility certificates if max( norm(dat.A * pt.x, Inf), - norm((pt.x - pt.xl) .* dat.lflag, Inf), - norm((pt.x + pt.xu) .* dat.uflag, Inf) + norm((pt.x .- pt.xl) .* dat.lflag, Inf), + norm((pt.x .+ pt.xu) .* dat.uflag, Inf) ) * (norm(dat.c, Inf) / max(1, norm(dat.b, Inf))) < - ϵi * dot(dat.c, pt.x) # Dual infeasible, i.e., primal unbounded mpc.primal_status = Sln_InfeasibilityCertificate @@ -197,7 +195,7 @@ function update_solver_status!(mpc::MPC{T}, ϵp::T, ϵd::T, ϵg::T, ϵi::T) wher return nothing end - δ = dat.A' * pt.y + (pt.zl .* dat.lflag) - (pt.zu .* dat.uflag) + δ = dat.A' * pt.y .+ (pt.zl .* dat.lflag) .- (pt.zu .* dat.uflag) if norm(δ, Inf) * max( norm(dat.l .* dat.lflag, Inf), norm(dat.u .* dat.uflag, Inf), @@ -367,11 +365,11 @@ function compute_starting_point(mpc::MPC{T}) where{T} # I. Recover positive primal-dual coordinates δx = one(T) + max( zero(T), - (-3 // 2) * minimum((pt.x - dat.l) .* dat.lflag), - (-3 // 2) * minimum((dat.u - pt.x) .* dat.uflag) + (-3 // 2) * minimum((pt.x .- dat.l) .* dat.lflag), + (-3 // 2) * minimum((dat.u .- pt.x) .* dat.uflag) ) - pt.xl .= ((pt.x - dat.l) .+ δx) .* dat.lflag - pt.xu .= ((dat.u - pt.x) .+ δx) .* dat.uflag + @. pt.xl = ((pt.x - dat.l) + δx) * dat.lflag + @. pt.xu = ((dat.u - pt.x) + δx) * dat.uflag z = dat.c - dat.A' * pt.y #= @@ -385,8 +383,8 @@ function compute_starting_point(mpc::MPC{T}) where{T} no | no | 0 | 0 | ----+-----+--------+---------+ =# - pt.zl .= ( z ./ (dat.lflag + dat.uflag)) .* dat.lflag - pt.zu .= (-z ./ (dat.lflag + dat.uflag)) .* dat.uflag + @. pt.zl = ( z / (dat.lflag + dat.uflag)) * dat.lflag + @. pt.zu = (-z / (dat.lflag + dat.uflag)) * dat.uflag δz = one(T) + max(zero(T), (-3 // 2) * minimum(pt.zl), (-3 // 2) * minimum(pt.zu)) pt.zl[dat.lflag] .+= δz diff --git a/src/IPM/MPC/step.jl b/src/IPM/MPC/step.jl index f1440aea..00f4264b 100644 --- a/src/IPM/MPC/step.jl +++ b/src/IPM/MPC/step.jl @@ -179,12 +179,10 @@ function solve_newton_system!(Δ::Point{T, Tv}, # II. Recover Δxl, Δxu @timeit mpc.timer "Δxl" begin - @. Δ.xl = -ξl + Δ.x - Δ.xl .*= dat.lflag + @. Δ.xl = (-ξl + Δ.x) * dat.lflag end @timeit mpc.timer "Δxu" begin - @. Δ.xu = ξu - Δ.x - Δ.xu .*= dat.uflag + @. Δ.xu = ( ξu - Δ.x) * dat.uflag end # III. Recover Δzl, Δzu @@ -259,8 +257,8 @@ function compute_corrector!(mpc::MPC{T, Tv}) where{T, Tv<:AbstractVector{T}} # Step length for affine-scaling direction αp_aff, αd_aff = mpc.αp, mpc.αd μₐ = ( - dot((pt.xl + αp_aff .* Δ.xl) .* dat.lflag, pt.zl + αd_aff .* Δ.zl) - + dot((pt.xu + αp_aff .* Δ.xu) .* dat.uflag, pt.zu + αd_aff .* Δ.zu) + dot((@. ((pt.xl + αp_aff * Δ.xl) * dat.lflag)), pt.zl .+ αd_aff .* Δ.zl) + + dot((@. ((pt.xu + αp_aff * Δ.xu) * dat.uflag)), pt.zu .+ αd_aff .* Δ.zu) ) / pt.p σ = clamp((μₐ / pt.μ)^3, sqrt(eps(T)), one(T) - sqrt(eps(T))) @@ -296,8 +294,8 @@ function compute_extra_correction!(mpc::MPC{T, Tv}; αd_ = min(αd + δ, one(T)) g = dot(pt.xl, pt.zl) + dot(pt.xu, pt.zu) - gₐ = dot((pt.xl + mpc.αp .* Δ.xl) .* dat.lflag, pt.zl + mpc.αd .* Δ.zl) + - dot((pt.xu + mpc.αp .* Δ.xu) .* dat.uflag, pt.zu + mpc.αd .* Δ.zu) + gₐ = dot((@. ((pt.xl + mpc.αp * Δ.xl) * dat.lflag)), pt.zl .+ mpc.αd .* Δ.zl) + + dot((@. ((pt.xu + mpc.αp * Δ.xu) * dat.uflag)), pt.zu .+ mpc.αd .* Δ.zu) μ = (gₐ / g) * (gₐ / g) * (gₐ / pt.p) # Newton RHS diff --git a/src/IPM/ipmdata.jl b/src/IPM/ipmdata.jl index 94999778..7222ada0 100644 --- a/src/IPM/ipmdata.jl +++ b/src/IPM/ipmdata.jl @@ -130,7 +130,7 @@ function IPMData(pb::ProblemData{T}, mfact::Factory) where{T} c0 = pb.obj0 if !pb.objsense # Flip objective for maximization problem - c .= -c + c .= .-c c0 = -c0 end @@ -170,4 +170,4 @@ function IPMData(pb::ProblemData{T}, mfact::Factory) where{T} u = [pb.uvar; uslack] return IPMData(A, b, pb.objsense, c, c0, l, u) -end \ No newline at end of file +end