[perf] Fuse more vectorized operations (#127)

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,16 +171,16 @@ 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
         hsd.solver_status = Trm_DualInfeasible
         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),

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

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,16 +186,16 @@ 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
         mpc.solver_status = Trm_DualInfeasible
         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

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

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
+end