From ac4b6458a4f6b5157f0fd4699af0db986ab70342 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 13:34:22 +1300
Subject: [PATCH 1/9] [docs] add a tutorial for PDHG and a basic MOI interface

---
 docs/make.jl                          |   1 +
 docs/src/tutorials/algorithms/pdhg.jl | 329 ++++++++++++++++++++++++++
 2 files changed, 330 insertions(+)
 create mode 100644 docs/src/tutorials/algorithms/pdhg.jl

diff --git a/docs/make.jl b/docs/make.jl
index be4b2f069cb..a3afbdfee9b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -380,6 +380,7 @@ const _PAGES = [
             "tutorials/algorithms/tsp_lazy_constraints.md",
             "tutorials/algorithms/rolling_horizon.md",
             "tutorials/algorithms/parallelism.md",
+            "tutorials/algorithms/pdhg.md",
         ],
         "Applications" => [
             "tutorials/applications/power_systems.md",
diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
new file mode 100644
index 00000000000..cefede83865
--- /dev/null
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -0,0 +1,329 @@
+# Copyright (c) 2024 Oscar Dowson and contributors                               #src
+#                                                                                #src
+# Permission is hereby granted, free of charge, to any person obtaining a copy   #src
+# of this software and associated documentation files (the "Software"), to deal  #src
+# in the Software without restriction, including without limitation the rights   #src
+# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell      #src
+# copies of the Software, and to permit persons to whom the Software is          #src
+# furnished to do so, subject to the following conditions:                       #src
+#                                                                                #src
+# The above copyright notice and this permission notice shall be included in all #src
+# copies or substantial portions of the Software.                                #src
+#                                                                                #src
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR     #src
+# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,       #src
+# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE    #src
+# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER         #src
+# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,  #src
+# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE  #src
+# SOFTWARE.                                                                      #src
+
+# # Writing a solver interface
+
+# The purpose of this tutorial is to demonstrate
+
+using JuMP
+import LinearAlgebra
+import MathOptInterface as MOI
+import Printf
+import SparseArrays
+
+# ## PDHG
+
+"""
+    solve_pdhg(
+        A::SparseArrays.SparseMatrixCSC{Float64,Int},
+        b::Vector{Float64},
+        c::Vector{Float64};
+        maximum_iterations::Int = 10_000,
+        tol::Float64 = 1e-4,
+        verbose::Bool = true,
+        log_frequency::Int = 1_000,
+    ) -> status, iterations, x, y
+
+A pedagogical implementation of PDHG that solves the linear program:
+
+```math
+\\begin{aligned}
+\\text{min}        \\ & c^\\top x \\
+\\text{subject to} \\ & Ax = b \\
+                      & x \\ge 0.
+\\end{aligned}
+```
+
+Note that this implementation is intentionally kept simple. It is not robust nor
+efficient, and it does not incorporate the theoretical improvements in the PDLP
+paper.
+
+## Keyword arguments
+
+ * `maximum_iterations::Int = 10_000`: the maximum number of iterations before
+   termination.
+ * `tol::Float64 = 1e-4`: the combined primal, dual, and strong duality
+   tolerance.
+ * `verbose::Bool = true`: print iteration logs
+ * `log_frequency::Int = 1_000`: print iteration logs every N iterations
+"""
+function solve_pdhg(
+    A::SparseArrays.SparseMatrixCSC{Float64,Int},
+    b::Vector{Float64},
+    c::Vector{Float64};
+    maximum_iterations::Int = 100_000,
+    tol::Float64 = 1e-4,
+    verbose::Bool = true,
+    log_frequency::Int = 1_000,
+)
+    printf(x::Float64) = Printf.@sprintf("% 1.6e", x)
+    printf(x::Int) = Printf.@sprintf("%6d", x)
+    m, n = size(A)
+    η = τ = 1 / LinearAlgebra.norm(A) - 1e-6
+    x, y, k, status = zeros(n), zeros(m), 0, MOI.OTHER_ERROR
+    if verbose
+        println(
+            "  iter      pobj          dobj         pfeas         dfeas       objfeas",
+        )
+    end
+    while status == MOI.OTHER_ERROR
+        x_next = max.(0.0, x - η * (A' * y + c))
+        y += τ * (A * (2 * x_next - x) - b)
+        x = x_next
+        k += 1
+        pfeas = LinearAlgebra.norm(A * x - b)
+        dfeas = LinearAlgebra.norm(max.(0.0, -A' * y - c))
+        objfeas = abs(c' * x + b' * y)
+        if pfeas <= tol && dfeas <= tol && objfeas <= tol
+            status = MOI.OPTIMAL
+        elseif k == maximum_iterations
+            status = MOI.ITERATION_LIMIT
+        end
+        if verbose && (mod(k, log_frequency) == 0 || status != MOI.OTHER_ERROR)
+            logs = printf.((k, c' * x, -b' * y, pfeas, dfeas, objfeas))
+            println(join(logs, " "))
+        end
+    end
+    return status, k, x, y
+end
+
+# Here's an example:
+
+A = [0.0 -1.0 -1.0 0.0 0.0; 6.0 8.0 0.0 -1.0 0.0; 7.0 12.0 0.0 0.0 -1.0]
+b = [-3.0, 100.0, 120.0]
+c = [12.0, 20.0, 0.0, 0.0, 0.0]
+status, k, x, y = solve_pdhg(SparseArrays.sparse(A), b, c)
+
+# ## The MOI interface
+
+# Converting example linear program from the modeler's form into the standard
+# form required by PDHG is tedious and error-prone. This section walks through
+# how to implement a basic interface to MathOptInterface, so that we can use our
+# algorithm from JuMP.
+
+"""
+    Optimizer()
+
+Create a new optimizer for PDHG.
+"""
+mutable struct Optimizer <: MOI.AbstractOptimizer
+    ## Information from solve_pdhg
+    status::MOI.TerminationStatusCode
+    iterations::Int
+    x::Vector{Float64}
+    y::Vector{Float64}
+    ## Other useful quantities
+    solve_time::Float64
+    obj_value::Float64
+
+    function Optimizer()
+        return new(MOI.OPTIMIZE_NOT_CALLED, 0, Float64[], Float64[], 0.0, 0.0)
+    end
+end
+
+# First, we need to implement two methods: [`MOI.is_empty`](@ref) and
+# [`MOI.empty!`](@ref). These are called whenever MOI needs to ensure that the
+# optimizer is in a clean state.
+
+function MOI.is_empty(model::Optimizer)
+    ## You might want to check every field, not just the status
+    return model.status == MOI.OPTIMIZE_NOT_CALLED
+end
+
+function MOI.empty!(model::Optimizer)
+    model.status = MOI.OPTIMIZE_NOT_CALLED
+    model.iterations = 0
+    model.solve_time = 0.0
+    model.obj_value = 0.0
+    empty!(model.x)
+    empty!(model.y)
+    return
+end
+
+# Next, we need to define what constraints the optimizer supports. Since our
+# standard form was $Ax = b$, we support only $Ax + b \in \{0\}$, which is a
+# [`MOI.VectorAffineFunction`](@ref) in [`MOI.Zeros`](@ref) constraint. Note
+# that you might have expected $Ax - b \in \{0\}$. We'll address the difference
+# in the sign of `b` in a few places later on.
+
+function MOI.supports_constraint(
+    ::Optimizer,
+    ::Type{MOI.VectorAffineFunction{Float64}},
+    ::Type{MOI.Zeros},
+)
+    return true
+end
+
+# By default, MOI assumes that it can add free variables. This isn't true for
+# our standard form, because we support only $x \ge 0$. Let's tell MOI that:
+
+MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false
+
+function MOI.supports_add_constrained_variables(
+    ::Optimizer,
+    ::Type{MOI.Nonnegatives},
+)
+    return true
+end
+
+# Finally, the objective function that we support is
+# [`MOI.ScalarAffineFunction`](@ref):
+
+function MOI.supports(
+    ::Optimizer,
+    ::MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}},
+)
+    return true
+end
+
+# Finally, we'll implement [`MOI.SolverName`](@ref) so that MOI knows how to
+# print the name of our optimizer:
+
+MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"
+
+
+MOI.Utilities.@product_of_sets(LinearZero, MOI.Zeros)
+
+const Cache = MOI.Utilities.GenericModel{
+    ## The coefficient type is Float64
+    Float64,
+    ## We use the default objective container
+    MOI.Utilities.ObjectiveContainer{Float64},
+    ## We use the default variable container
+    MOI.Utilities.VariablesContainer{Float64},
+    ## We use a Matrix of Constraints to represent `A * x + b in K`
+    MOI.Utilities.MatrixOfConstraints{
+        ## The number type is Float64
+        Float64,
+        ## The matrix type `A` is a sparse matrix
+        MOI.Utilities.MutableSparseMatrixCSC{
+            ## ... with Float64 coefficients
+            Float64,
+            ## ... Int64 row and column indices
+            Int,
+            ## ... and it uses one-based indexing
+            MOI.Utilities.OneBasedIndexing,
+        },
+        ## The vector type `b` is a Julia `Vector`
+        Vector{Float64},
+        ## The set type `K` is the LinearZero set we defined above
+        LinearZero{Float64},
+    },
+}
+
+# If you were interfacing with a solver written in C that expected zero-based
+# indices, you might use instead:
+MOI.Utilities.MutableSparseMatrixCSC{
+    Cdouble,
+    Cint,
+    MOI.Utilities.ZeroBasedIndexing,
+}
+
+function MOI.add_constrained_variables(model::Cache, set::MOI.Nonnegatives)
+    x = MOI.add_variables(model, MOI.dimension(set))
+    MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
+    ci = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}(x[1].value)
+    return x, ci
+end
+
+function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
+    start_time = time()
+    cache = Cache()
+    index_map = MOI.copy_to(cache, src)
+    A = convert(
+        SparseArrays.SparseMatrixCSC{Float64,Int},
+        cache.constraints.coefficients,
+    )
+    ## MOI models Ax = b as Ax + b in {0}, so b differs by -
+    b = -cache.constraints.constants
+    c = zeros(size(A, 2))
+    sense = ifelse(cache.objective.sense == MOI.MAX_SENSE, -1, 1)
+    F = MOI.ScalarAffineFunction{Float64}
+    obj = MOI.get(src, MOI.ObjectiveFunction{F}())
+    for term in obj.terms
+        c[term.variable.value] += sense * term.coefficient
+    end
+    dest.status, dest.iterations, dest.x, dest.y = solve_pdhg(A, b, c)
+    dest.obj_value = obj.constant + c' * dest.x
+    dest.solve_time = time() - start_time
+    return index_map, false
+end
+
+
+function MOI.get(model::Optimizer, ::MOI.ResultCount)
+    return model.status == MOI.OPTIMAL ? 1 : 0
+end
+
+MOI.get(model::Optimizer, ::MOI.RawStatusString) = string(model.status)
+
+MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.status
+
+function MOI.get(model::Optimizer, attr::Union{MOI.PrimalStatus,MOI.DualStatus})
+    if attr.result_index == 1 && model.status == MOI.OPTIMAL
+        return MOI.FEASIBLE_POINT
+    end
+    return MOI.NO_SOLUTION
+end
+
+# Now we can implement [`MOI.ObjectiveValue`](@ref), [`MOI.VariablePrimal`](@ref),
+# and [`MOI.ConstraintDual`](@ref):
+
+function MOI.get(model::Optimizer, attr::MOI.ObjectiveValue)
+    MOI.check_result_index_bounds(model, attr)
+    return model.obj_value
+end
+
+function MOI.get(
+    model::Optimizer,
+    attr::MOI.VariablePrimal,
+    x::MOI.VariableIndex,
+)
+    MOI.check_result_index_bounds(model, attr)
+    return model.x[x.value]
+end
+
+function MOI.get(
+    model::Optimizer,
+    attr::MOI.ConstraintDual,
+    ci::MOI.ConstraintIndex,
+)
+    MOI.check_result_index_bounds(model, attr)
+    ## MOI models Ax = b as Ax + b in {0}, so the dual differs by -
+    return -model.y[ci.value]
+end
+
+# Some other useful result quantities are [`MOI.SolveTimeSec`](@ref) and
+# [`MOI.BarrierIterations`](@ref):
+
+MOI.get(model::Optimizer, ::MOI.SolveTimeSec) = model.solve_time
+MOI.get(model::Optimizer, ::MOI.BarrierIterations) = model.iterations
+
+# ## A JuMP example
+
+# Now we can solve an arbitrary linear program with JuMP:
+
+model = Model(Optimizer)
+@variable(model, x >= 0)
+@variable(model, 0 <= y <= 3)
+@objective(model, Min, 12x + 20y)
+@constraint(model, c1, 6x + 8y >= 100)
+@constraint(model, c2, 7x + 12y >= 120)
+optimize!(model)
+solution_summary(model; verbose = true)

From 2c43bc8cdff2ad246b62ef7273a6de3eeaba27cb Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 14:15:42 +1300
Subject: [PATCH 2/9] Update

---
 docs/src/tutorials/algorithms/pdhg.jl | 2 --
 1 file changed, 2 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index cefede83865..ee6a5fd42a2 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -198,7 +198,6 @@ end
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"
 
-
 MOI.Utilities.@product_of_sets(LinearZero, MOI.Zeros)
 
 const Cache = MOI.Utilities.GenericModel{
@@ -266,7 +265,6 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     return index_map, false
 end
 
-
 function MOI.get(model::Optimizer, ::MOI.ResultCount)
     return model.status == MOI.OPTIMAL ? 1 : 0
 end

From 2d1ee02e7583fa39ea5c57160175578ae2b60809 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 15:09:24 +1300
Subject: [PATCH 3/9] Update

---
 docs/src/tutorials/algorithms/pdhg.jl | 160 +++++++++++++++++---------
 1 file changed, 106 insertions(+), 54 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index ee6a5fd42a2..be12223c452 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -20,7 +20,11 @@
 
 # # Writing a solver interface
 
-# The purpose of this tutorial is to demonstrate
+# The purpose of this tutorial is to demonstrate how to implement a basic solver
+# interface to MathOptInterface. As a motivating example, we implement the
+# Primal Dual Hybrid Gradient (PDHG) method. PDHG is a first-order method that
+# can solve convex optimization problems.
+# Google has a [good introduction to the math behind PDLP](https://developers.google.com/optimization/lp/pdlp_math).
 
 using JuMP
 import LinearAlgebra
@@ -28,42 +32,21 @@ import MathOptInterface as MOI
 import Printf
 import SparseArrays
 
-# ## PDHG
+# ## Primal Dual Hybrid Gradient
+
+# Here is a pedagogical implementation of PDHG that solves the linear program:
+# ```math
+# \\begin{aligned}
+# \\text{min}        \\ & c^\\top x \\
+# \\text{subject to} \\ & Ax = b \\
+#                       & x \\ge 0.
+# \\end{aligned}
+# ```
+#
+# Note that this implementation is intentionally kept simple. It is not robust
+# nor efficient, and it does not incorporate the theoretical improvements in the
+# PDLP paper.
 
-"""
-    solve_pdhg(
-        A::SparseArrays.SparseMatrixCSC{Float64,Int},
-        b::Vector{Float64},
-        c::Vector{Float64};
-        maximum_iterations::Int = 10_000,
-        tol::Float64 = 1e-4,
-        verbose::Bool = true,
-        log_frequency::Int = 1_000,
-    ) -> status, iterations, x, y
-
-A pedagogical implementation of PDHG that solves the linear program:
-
-```math
-\\begin{aligned}
-\\text{min}        \\ & c^\\top x \\
-\\text{subject to} \\ & Ax = b \\
-                      & x \\ge 0.
-\\end{aligned}
-```
-
-Note that this implementation is intentionally kept simple. It is not robust nor
-efficient, and it does not incorporate the theoretical improvements in the PDLP
-paper.
-
-## Keyword arguments
-
- * `maximum_iterations::Int = 10_000`: the maximum number of iterations before
-   termination.
- * `tol::Float64 = 1e-4`: the combined primal, dual, and strong duality
-   tolerance.
- * `verbose::Bool = true`: print iteration logs
- * `log_frequency::Int = 1_000`: print iteration logs every N iterations
-"""
 function solve_pdhg(
     A::SparseArrays.SparseMatrixCSC{Float64,Int},
     b::Vector{Float64},
@@ -109,14 +92,39 @@ end
 A = [0.0 -1.0 -1.0 0.0 0.0; 6.0 8.0 0.0 -1.0 0.0; 7.0 12.0 0.0 0.0 -1.0]
 b = [-3.0, 100.0, 120.0]
 c = [12.0, 20.0, 0.0, 0.0, 0.0]
-status, k, x, y = solve_pdhg(SparseArrays.sparse(A), b, c)
+status, k, x, y = solve_pdhg(SparseArrays.sparse(A), b, c);
+
+# The termination status is:
+
+status
+
+# The solve took the following number of iterations:
+
+k
+
+# The primal solution is:
+
+x
+
+# The dual multipliers are:
+
+y
 
 # ## The MOI interface
 
-# Converting example linear program from the modeler's form into the standard
-# form required by PDHG is tedious and error-prone. This section walks through
-# how to implement a basic interface to MathOptInterface, so that we can use our
-# algorithm from JuMP.
+# Converting example linear program from the modeler's form into the `A`, `b`,
+# and `c` matrices of the standard form required by our implementation of PDHG
+# is tedious and error-prone. This section walks through how to implement a
+# basic interface to MathOptInterface, so that we can use our algorithm from
+# JuMP.
+
+# ### The Optimizer type
+
+# Create an optimizer by subtyping [`MOI.AbstractOptimizer`](@ref). By
+# convention, the name of this type is `Optimizer`, and most optimizers are
+# available as `PackageName.Optimizer`.
+
+# The fields inside the optimizer are arbitrary. Store whatever is useful.
 
 """
     Optimizer()
@@ -138,9 +146,9 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     end
 end
 
-# First, we need to implement two methods: [`MOI.is_empty`](@ref) and
-# [`MOI.empty!`](@ref). These are called whenever MOI needs to ensure that the
-# optimizer is in a clean state.
+# Now that we have an `Optimizer`, we need to implement two methods:
+# [`MOI.is_empty`](@ref) and [`MOI.empty!`](@ref). These are called whenever MOI
+# needs to ensure that the optimizer is in a clean state.
 
 function MOI.is_empty(model::Optimizer)
     ## You might want to check every field, not just the status
@@ -198,9 +206,24 @@ end
 
 MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"
 
-MOI.Utilities.@product_of_sets(LinearZero, MOI.Zeros)
+# ### GenericModel
+
+# The simplest way to solve a problem with your optimizer is to implement the
+# method `MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)`, where `src` is an
+# input model and `dest` is your empty optimizer.
+
+# To implement this method you would need to query the variables and constraints
+# in `src` and the convert these into the matrix data expected by `solve_pdhg`.
+# Since matrix input is a common requirement of solvers, MOI includes utilities
+# to simplify the process.
 
-const Cache = MOI.Utilities.GenericModel{
+# The downside of the utility is that involves a highly parameterized type with
+# a large number of possible configurations.The upside of the utility is that,
+# once setup, it requires few lines of code to extract the constraint matrices.
+
+MOI.Utilities.@product_of_sets(SetOfZeros, MOI.Zeros)
+
+const CacheModel = MOI.Utilities.GenericModel{
     ## The coefficient type is Float64
     Float64,
     ## We use the default objective container
@@ -222,20 +245,26 @@ const Cache = MOI.Utilities.GenericModel{
         },
         ## The vector type `b` is a Julia `Vector`
         Vector{Float64},
-        ## The set type `K` is the LinearZero set we defined above
-        LinearZero{Float64},
+        ## The set type `K` is the SetOfZeros type we defined above
+        SetOfZeros{Float64},
     },
 }
 
-# If you were interfacing with a solver written in C that expected zero-based
-# indices, you might use instead:
+# As one example of possible alternate configuration, if you were interfacing
+# with a solver written in C that expected zero-based indices, you might use
+# instead:
+
 MOI.Utilities.MutableSparseMatrixCSC{
     Cdouble,
     Cint,
     MOI.Utilities.ZeroBasedIndexing,
 }
 
-function MOI.add_constrained_variables(model::Cache, set::MOI.Nonnegatives)
+# !!! tip
+#     The best place to look at how to configure `GenericModel` is to find an
+#     existing solver with the same input standard form that you require.
+
+function MOI.add_constrained_variables(model::CacheModel, set::MOI.Nonnegatives)
     x = MOI.add_variables(model, MOI.dimension(set))
     MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
     ci = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}(x[1].value)
@@ -244,7 +273,7 @@ end
 
 function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     start_time = time()
-    cache = Cache()
+    cache = CacheModel()
     index_map = MOI.copy_to(cache, src)
     A = convert(
         SparseArrays.SparseMatrixCSC{Float64,Int},
@@ -252,10 +281,10 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     )
     ## MOI models Ax = b as Ax + b in {0}, so b differs by -
     b = -cache.constraints.constants
-    c = zeros(size(A, 2))
     sense = ifelse(cache.objective.sense == MOI.MAX_SENSE, -1, 1)
     F = MOI.ScalarAffineFunction{Float64}
     obj = MOI.get(src, MOI.ObjectiveFunction{F}())
+    c = zeros(size(A, 2))
     for term in obj.terms
         c[term.variable.value] += sense * term.coefficient
     end
@@ -265,11 +294,31 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     return index_map, false
 end
 
+# ## Solutions
+
+# Now that we've solved a model, let's implement the required solution
+# attributes.
+
+# First, we need to tell MOI how many solutions we found via
+# [`MOI.ResultCount`](@ref):
+
 function MOI.get(model::Optimizer, ::MOI.ResultCount)
     return model.status == MOI.OPTIMAL ? 1 : 0
 end
 
-MOI.get(model::Optimizer, ::MOI.RawStatusString) = string(model.status)
+# and implement [`MOI.RawStatusString`](@ref) to provide a user-readable string
+# that describes what happened:
+
+function MOI.get(model::Optimizer, ::MOI.RawStatusString)
+    if model.status == MOI.OPTIMAL
+        return "found a primal-dual optimal solution (subject to tolerances)"
+    end
+    return "failed to solve"
+end
+
+# Then, we need to implement the three types of problem status:
+# [`MOI.TerminationStatus`](@ref), [`MOI.PrimalStatus`](@ref) and
+# [`MOI.DualStatus`](@ref):
 
 MOI.get(model::Optimizer, ::MOI.TerminationStatus) = model.status
 
@@ -324,4 +373,7 @@ model = Model(Optimizer)
 @constraint(model, c1, 6x + 8y >= 100)
 @constraint(model, c2, 7x + 12y >= 120)
 optimize!(model)
+
+#-
+
 solution_summary(model; verbose = true)

From a7e970f0d11663cc6a04c207f2d133cb3cca0cc7 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 15:39:35 +1300
Subject: [PATCH 4/9] Update

---
 docs/src/tutorials/algorithms/pdhg.jl | 88 ++++++++++++++++++++++++---
 1 file changed, 80 insertions(+), 8 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index be12223c452..6083fc4f9dd 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -132,6 +132,11 @@ y
 Create a new optimizer for PDHG.
 """
 mutable struct Optimizer <: MOI.AbstractOptimizer
+    x_to_col::Dict{MOI.VariableIndex,Int}
+    ci_to_rows::Dict{
+        MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},MOI.Zeros},
+        Vector{Int},
+    }
     ## Information from solve_pdhg
     status::MOI.TerminationStatusCode
     iterations::Int
@@ -142,7 +147,17 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     obj_value::Float64
 
     function Optimizer()
-        return new(MOI.OPTIMIZE_NOT_CALLED, 0, Float64[], Float64[], 0.0, 0.0)
+        F = MOI.VectorAffineFunction{Float64}
+        return new(
+            Dict{MOI.VariableIndex,Int}(),
+            Dict{MOI.ConstraintIndex{F,MOI.Zeros},Vector{Int}}(),
+            MOI.OPTIMIZE_NOT_CALLED,
+            0,
+            Float64[],
+            Float64[],
+            0.0,
+            0.0,
+        )
     end
 end
 
@@ -151,11 +166,13 @@ end
 # needs to ensure that the optimizer is in a clean state.
 
 function MOI.is_empty(model::Optimizer)
-    ## You might want to check every field, not just the status
-    return model.status == MOI.OPTIMIZE_NOT_CALLED
+    ## You might want to check every field, not just a few
+    return isempty(model.x_to_col) && model.status == MOI.OPTIMIZE_NOT_CALLED
 end
 
 function MOI.empty!(model::Optimizer)
+    empty!(model.x_to_col)
+    empty!(model.ci_to_rows)
     model.status = MOI.OPTIMIZE_NOT_CALLED
     model.iterations = 0
     model.solve_time = 0.0
@@ -221,6 +238,8 @@ MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"
 # a large number of possible configurations.The upside of the utility is that,
 # once setup, it requires few lines of code to extract the constraint matrices.
 
+# Define the set of sets that our standard form supports. For PDHG, we support
+# only `Ax + b in {0}`:
 MOI.Utilities.@product_of_sets(SetOfZeros, MOI.Zeros)
 
 const CacheModel = MOI.Utilities.GenericModel{
@@ -264,6 +283,10 @@ MOI.Utilities.MutableSparseMatrixCSC{
 #     The best place to look at how to configure `GenericModel` is to find an
 #     existing solver with the same input standard form that you require.
 
+# We need to make one modification to `CacheModel` to tell MOI that
+# $x \in \mathbb{R}_+$ is equivalent to adding variables in
+# [`MOI.GreaterThan`](@ref):
+
 function MOI.add_constrained_variables(model::CacheModel, set::MOI.Nonnegatives)
     x = MOI.add_variables(model, MOI.dimension(set))
     MOI.add_constraint.(model, x, MOI.GreaterThan(0.0))
@@ -271,16 +294,32 @@ function MOI.add_constrained_variables(model::CacheModel, set::MOI.Nonnegatives)
     return x, ci
 end
 
+# ### The optimize method
+
 function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     start_time = time()
     cache = CacheModel()
     index_map = MOI.copy_to(cache, src)
+    ## Create a map of the constraint indices to their row in the `dest` matrix
+    F, S = MOI.VectorAffineFunction{Float64}, MOI.Zeros
+    for src_ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
+        dest_ci = index_map[src_ci]
+        dest.ci_to_rows[dest_ci] =
+            MOI.Utilities.rows(cache.constraints.sets, dest_ci)
+    end
+    ## Create a map of the variable indices to their column in the `dest` matrix
+    for (i, src_x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
+        dest.x_to_col[index_map[src_x]] = i
+    end
+    ## Now we can access the `A` matrix:
     A = convert(
         SparseArrays.SparseMatrixCSC{Float64,Int},
         cache.constraints.coefficients,
     )
     ## MOI models Ax = b as Ax + b in {0}, so b differs by -
     b = -cache.constraints.constants
+    ## The `c` vector is more involved, because we need to account for the
+    ## objective sense:
     sense = ifelse(cache.objective.sense == MOI.MAX_SENSE, -1, 1)
     F = MOI.ScalarAffineFunction{Float64}
     obj = MOI.get(src, MOI.ObjectiveFunction{F}())
@@ -288,9 +327,14 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     for term in obj.terms
         c[term.variable.value] += sense * term.coefficient
     end
+    ## Now we can solve the problem with PDHG and record the solution:
     dest.status, dest.iterations, dest.x, dest.y = solve_pdhg(A, b, c)
-    dest.obj_value = obj.constant + c' * dest.x
+    ## as well as record two derived quantities: the objective value and the
+    ## solve time
+    dest.obj_value = obj.constant + sense * c' * dest.x
     dest.solve_time = time() - start_time
+    ## We need to return the index map, and `false`, to indicate to MOI that we
+    ## do not support incremental modification of the model.
     return index_map, false
 end
 
@@ -343,17 +387,17 @@ function MOI.get(
     x::MOI.VariableIndex,
 )
     MOI.check_result_index_bounds(model, attr)
-    return model.x[x.value]
+    return model.x[model.x_to_col[x]]
 end
 
 function MOI.get(
     model::Optimizer,
     attr::MOI.ConstraintDual,
-    ci::MOI.ConstraintIndex,
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},MOI.Zeros},
 )
     MOI.check_result_index_bounds(model, attr)
     ## MOI models Ax = b as Ax + b in {0}, so the dual differs by -
-    return -model.y[ci.value]
+    return -model.y[model.ci_to_rows[ci]]
 end
 
 # Some other useful result quantities are [`MOI.SolveTimeSec`](@ref) and
@@ -364,7 +408,20 @@ MOI.get(model::Optimizer, ::MOI.BarrierIterations) = model.iterations
 
 # ## A JuMP example
 
-# Now we can solve an arbitrary linear program with JuMP:
+# Now we can solve an arbitrary linear program with JuMP. Here's the same
+# standard form as before:
+
+model = Model(Optimizer)
+@variable(model, x[1:5] >= 0)
+@objective(model, Min, c' * x)
+@constraint(model, c3, A * x == b)
+optimize!(model)
+
+#-
+
+solution_summary(model; verbose = true)
+
+# But we could also have written:
 
 model = Model(Optimizer)
 @variable(model, x >= 0)
@@ -377,3 +434,18 @@ optimize!(model)
 #-
 
 solution_summary(model; verbose = true)
+
+# Other variations are also possible:
+
+model = Model(Optimizer)
+@variable(model, x[1:5] >= 0)
+@objective(model, Max, -c' * x)
+@constraint(model, c4, A * x .== b)
+optimize!(model)
+
+#-
+
+solution_summary(model; verbose = true)
+
+# Behind the scenes, JuMP and MathOptInterface reformulate the problem from the
+# modeller's form into the standard form defined by our `Optimizer`.

From 2dab539cf5adc87e9aa30fd9e3d19933f8bd2857 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 15:52:30 +1300
Subject: [PATCH 5/9] Fix latex

---
 docs/src/tutorials/algorithms/pdhg.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index 6083fc4f9dd..b162c416ccf 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -36,11 +36,11 @@ import SparseArrays
 
 # Here is a pedagogical implementation of PDHG that solves the linear program:
 # ```math
-# \\begin{aligned}
-# \\text{min}        \\ & c^\\top x \\
-# \\text{subject to} \\ & Ax = b \\
-#                       & x \\ge 0.
-# \\end{aligned}
+# \begin{aligned}
+# \text{min}        \ & c^\top x   \\
+# \text{subject to} \ & Ax = b     \\
+#                     & x \ge 0.
+# \end{aligned}
 # ```
 #
 # Note that this implementation is intentionally kept simple. It is not robust

From c1b2e31ae4c6487814b09811924854eb4c41d3b4 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Wed, 9 Oct 2024 17:04:57 +1300
Subject: [PATCH 6/9] Update

---
 docs/src/tutorials/algorithms/pdhg.jl | 86 +++++++++++++++++----------
 1 file changed, 55 insertions(+), 31 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index b162c416ccf..67c22ce0043 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -24,7 +24,13 @@
 # interface to MathOptInterface. As a motivating example, we implement the
 # Primal Dual Hybrid Gradient (PDHG) method. PDHG is a first-order method that
 # can solve convex optimization problems.
-# Google has a [good introduction to the math behind PDLP](https://developers.google.com/optimization/lp/pdlp_math).
+#
+# Google has a [good introduction to the math behind PDLP](https://developers.google.com/optimization/lp/pdlp_math),
+# which is a variant of PDHG specialized for linear programs.
+
+# ## Required packages
+
+# This tutorial requires the following packages:
 
 using JuMP
 import LinearAlgebra
@@ -34,7 +40,8 @@ import SparseArrays
 
 # ## Primal Dual Hybrid Gradient
 
-# Here is a pedagogical implementation of PDHG that solves the linear program:
+# The following function is a pedagogical implementation of PDHG that solves the
+# linear program:
 # ```math
 # \begin{aligned}
 # \text{min}        \ & c^\top x   \\
@@ -73,7 +80,7 @@ function solve_pdhg(
         k += 1
         pfeas = LinearAlgebra.norm(A * x - b)
         dfeas = LinearAlgebra.norm(max.(0.0, -A' * y - c))
-        objfeas = abs(c' * x + b' * y)
+        objfeas = abs(c' * x - -b' * y)
         if pfeas <= tol && dfeas <= tol && objfeas <= tol
             status = MOI.OPTIMAL
         elseif k == maximum_iterations
@@ -112,11 +119,12 @@ y
 
 # ## The MOI interface
 
-# Converting example linear program from the modeler's form into the `A`, `b`,
-# and `c` matrices of the standard form required by our implementation of PDHG
-# is tedious and error-prone. This section walks through how to implement a
-# basic interface to MathOptInterface, so that we can use our algorithm from
-# JuMP.
+# Converting a linear program from the modeler's form into the `A`, `b`, and `c`
+# matrices of the standard form required by our implementation of PDHG is
+# tedious and error-prone. This section walks through how to implement a basic
+# interface to MathOptInterface, so that we can use our algorithm from JuMP.
+
+# For a more comprehensive guide, see [Implementing a solver interface](@ref).
 
 # ### The Optimizer type
 
@@ -132,7 +140,9 @@ y
 Create a new optimizer for PDHG.
 """
 mutable struct Optimizer <: MOI.AbstractOptimizer
+    ## A mapping from variable to column
     x_to_col::Dict{MOI.VariableIndex,Int}
+    ## A mapping from constraint to rows
     ci_to_rows::Dict{
         MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64},MOI.Zeros},
         Vector{Int},
@@ -147,7 +157,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     obj_value::Float64
 
     function Optimizer()
-        F = MOI.VectorAffineFunction{Float64}
+        F =
         return new(
             Dict{MOI.VariableIndex,Int}(),
             Dict{MOI.ConstraintIndex{F,MOI.Zeros},Vector{Int}}(),
@@ -208,8 +218,7 @@ function MOI.supports_add_constrained_variables(
     return true
 end
 
-# Finally, the objective function that we support is
-# [`MOI.ScalarAffineFunction`](@ref):
+# The objective function that we support is [`MOI.ScalarAffineFunction`](@ref):
 
 function MOI.supports(
     ::Optimizer,
@@ -234,14 +243,19 @@ MOI.get(::Optimizer, ::MOI.SolverName) = "PDHG"
 # Since matrix input is a common requirement of solvers, MOI includes utilities
 # to simplify the process.
 
-# The downside of the utility is that involves a highly parameterized type with
-# a large number of possible configurations.The upside of the utility is that,
-# once setup, it requires few lines of code to extract the constraint matrices.
+# The downside of the utilities is that they involve a highly parameterized type
+# with a large number of possible configurations.The upside of the utilities is
+# that, once setup, they requires few lines of code to extract the problem
+# matrices.
+
+# First, we need to define the set of sets that our standard form supports. For
+# PDHG, we support only `Ax + b in {0}`:
 
-# Define the set of sets that our standard form supports. For PDHG, we support
-# only `Ax + b in {0}`:
 MOI.Utilities.@product_of_sets(SetOfZeros, MOI.Zeros)
 
+# Then, we define a [`MOI.Utilities.GenericModel`](@ref). This is the highly
+# parameterized type that can be customized.
+
 const CacheModel = MOI.Utilities.GenericModel{
     ## The coefficient type is Float64
     Float64,
@@ -296,27 +310,25 @@ end
 
 # ### The optimize method
 
+# Now we define the most important method for our optimizer.
+
 function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
+    ## In addition to the values returned by `solve_pdhg`, it may be useful to
+    ## record other attributes, such as the solve time.
     start_time = time()
+    ## Construct a cache to store our problem data:
     cache = CacheModel()
+    ## MOI includes a utility to copy an arbitrary `src` model into `cache`. The
+    ## return, `index_map`, is a mapping from indices in `src` to indices in
+    ## `dest`.
     index_map = MOI.copy_to(cache, src)
-    ## Create a map of the constraint indices to their row in the `dest` matrix
-    F, S = MOI.VectorAffineFunction{Float64}, MOI.Zeros
-    for src_ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
-        dest_ci = index_map[src_ci]
-        dest.ci_to_rows[dest_ci] =
-            MOI.Utilities.rows(cache.constraints.sets, dest_ci)
-    end
-    ## Create a map of the variable indices to their column in the `dest` matrix
-    for (i, src_x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
-        dest.x_to_col[index_map[src_x]] = i
-    end
     ## Now we can access the `A` matrix:
     A = convert(
         SparseArrays.SparseMatrixCSC{Float64,Int},
         cache.constraints.coefficients,
     )
-    ## MOI models Ax = b as Ax + b in {0}, so b differs by -
+    ## and the b vector (note that MOI models Ax = b as Ax + b in {0}, so b
+    ## differs by -):
     b = -cache.constraints.constants
     ## The `c` vector is more involved, because we need to account for the
     ## objective sense:
@@ -329,8 +341,20 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     end
     ## Now we can solve the problem with PDHG and record the solution:
     dest.status, dest.iterations, dest.x, dest.y = solve_pdhg(A, b, c)
-    ## as well as record two derived quantities: the objective value and the
-    ## solve time
+    ## To help assign the values of the x and y vectors to the appropriate
+    ## variables and constrats, we need a map of the constraint indices to their
+    ## row in the `dest` matrix and a map of the variable indices to their
+    # column in the `dest` matrix:
+    F, S = MOI.VectorAffineFunction{Float64}, MOI.Zeros
+    for src_ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
+        dest.ci_to_rows[index_map[src_ci]] =
+            MOI.Utilities.rows(cache.constraints.sets, index_map[src_ci])
+    end
+    for (i, src_x) in enumerate(MOI.get(src, MOI.ListOfVariableIndices()))
+        dest.x_to_col[index_map[src_x]] = i
+    end
+    ## We can now record two derived quantities: the primal objective value and
+    ## the solve time.
     dest.obj_value = obj.constant + sense * c' * dest.x
     dest.solve_time = time() - start_time
     ## We need to return the index map, and `false`, to indicate to MOI that we
@@ -340,7 +364,7 @@ end
 
 # ## Solutions
 
-# Now that we've solved a model, let's implement the required solution
+# Now that we know how to solve a model, let's implement the required solution
 # attributes.
 
 # First, we need to tell MOI how many solutions we found via

From 72c44bd9b6902f23ab3ef63a2ae8879c54668af7 Mon Sep 17 00:00:00 2001
From: Oscar Dowson <odow@users.noreply.github.com>
Date: Thu, 10 Oct 2024 08:31:40 +1300
Subject: [PATCH 7/9] Update docs/src/tutorials/algorithms/pdhg.jl

---
 docs/src/tutorials/algorithms/pdhg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index 67c22ce0043..0e2c3601e0b 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -157,7 +157,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     obj_value::Float64
 
     function Optimizer()
-        F =
+        F = MOI.VectorAffineFunction{Float64}
         return new(
             Dict{MOI.VariableIndex,Int}(),
             Dict{MOI.ConstraintIndex{F,MOI.Zeros},Vector{Int}}(),

From b050897dba20ddf02cfbdd34a4abe4025bb971ea Mon Sep 17 00:00:00 2001
From: Oscar Dowson <odow@users.noreply.github.com>
Date: Thu, 10 Oct 2024 09:23:21 +1300
Subject: [PATCH 8/9] Update docs/src/tutorials/algorithms/pdhg.jl

---
 docs/src/tutorials/algorithms/pdhg.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index 0e2c3601e0b..6302f2a6f8f 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -344,7 +344,7 @@ function MOI.optimize!(dest::Optimizer, src::MOI.ModelLike)
     ## To help assign the values of the x and y vectors to the appropriate
     ## variables and constrats, we need a map of the constraint indices to their
     ## row in the `dest` matrix and a map of the variable indices to their
-    # column in the `dest` matrix:
+    ## column in the `dest` matrix:
     F, S = MOI.VectorAffineFunction{Float64}, MOI.Zeros
     for src_ci in MOI.get(src, MOI.ListOfConstraintIndices{F,S}())
         dest.ci_to_rows[index_map[src_ci]] =

From 41d368e38f0e0d4c5cf516f65a27418d05a070f8 Mon Sep 17 00:00:00 2001
From: odow <o.dowson@gmail.com>
Date: Thu, 10 Oct 2024 10:24:02 +1300
Subject: [PATCH 9/9] Update

---
 docs/src/tutorials/algorithms/pdhg.jl | 39 +++++++++++++++++++++------
 1 file changed, 31 insertions(+), 8 deletions(-)

diff --git a/docs/src/tutorials/algorithms/pdhg.jl b/docs/src/tutorials/algorithms/pdhg.jl
index 6302f2a6f8f..818e795d64f 100644
--- a/docs/src/tutorials/algorithms/pdhg.jl
+++ b/docs/src/tutorials/algorithms/pdhg.jl
@@ -52,7 +52,8 @@ import SparseArrays
 #
 # Note that this implementation is intentionally kept simple. It is not robust
 # nor efficient, and it does not incorporate the theoretical improvements in the
-# PDLP paper.
+# PDLP paper. It does use two workspace vectors so that the body of the
+# iteration loop is non-allocating.
 
 function solve_pdhg(
     A::SparseArrays.SparseMatrixCSC{Float64,Int},
@@ -67,20 +68,42 @@ function solve_pdhg(
     printf(x::Int) = Printf.@sprintf("%6d", x)
     m, n = size(A)
     η = τ = 1 / LinearAlgebra.norm(A) - 1e-6
-    x, y, k, status = zeros(n), zeros(m), 0, MOI.OTHER_ERROR
+    x, x_next, y, k, status = zeros(n), zeros(n), zeros(m), 0, MOI.OTHER_ERROR
+    m_workspace, n_workspace = zeros(m), zeros(n)
     if verbose
         println(
             "  iter      pobj          dobj         pfeas         dfeas       objfeas",
         )
     end
     while status == MOI.OTHER_ERROR
-        x_next = max.(0.0, x - η * (A' * y + c))
-        y += τ * (A * (2 * x_next - x) - b)
-        x = x_next
         k += 1
-        pfeas = LinearAlgebra.norm(A * x - b)
-        dfeas = LinearAlgebra.norm(max.(0.0, -A' * y - c))
-        objfeas = abs(c' * x - -b' * y)
+        ## =====================================================================
+        ## This block computes x_next = max.(0.0, x - η * (A' * y + c))
+        LinearAlgebra.mul!(x_next, A', y)
+        LinearAlgebra.axpby!(-η, c, -η, x_next)
+        x_next .+= x
+        x_next .= max.(0.0, x_next)
+        ## =====================================================================
+        ## This block computes y += τ * (A * (2 * x_next - x) - b)
+        copy!(n_workspace, x_next)
+        LinearAlgebra.axpby!(-1.0, x, 2.0, n_workspace)
+        LinearAlgebra.mul!(y, A, n_workspace, τ, 1.0)
+        LinearAlgebra.axpy!(-τ, b, y)
+        ## =====================================================================
+        copy!(x, x_next)
+        ## =====================================================================
+        ## This block computes pfeas = LinearAlgebra.norm(A * x - b)
+        LinearAlgebra.mul!(m_workspace, A, x)
+        m_workspace .-= b
+        pfeas = LinearAlgebra.norm(m_workspace)
+        ## =====================================================================
+        ## This block computes dfeas = LinearAlgebra.norm(min.(0.0, A' * y + c))
+        LinearAlgebra.mul!(n_workspace, A', y)
+        n_workspace .+= c
+        n_workspace .= min.(0.0, n_workspace)
+        dfeas = LinearAlgebra.norm(n_workspace)
+        ## =====================================================================
+        objfeas = abs(LinearAlgebra.dot(c, x) + LinearAlgebra.dot(b, y))
         if pfeas <= tol && dfeas <= tol && objfeas <= tol
             status = MOI.OPTIMAL
         elseif k == maximum_iterations