diff --git a/docs/src/submodules/Test/reference.md b/docs/src/submodules/Test/reference.md
index 23a6740edd..830508f82f 100644
--- a/docs/src/submodules/Test/reference.md
+++ b/docs/src/submodules/Test/reference.md
@@ -18,4 +18,5 @@ Test.setup_test
 Test.version_added
 Test.@requires
 Test.RequirementUnmet
+Test.HS071
 ```
diff --git a/src/Test/test_nonlinear.jl b/src/Test/test_nonlinear.jl
index 575c2af642..81826b09dd 100644
--- a/src/Test/test_nonlinear.jl
+++ b/src/Test/test_nonlinear.jl
@@ -10,15 +10,18 @@
         enable_hessian_vector_product::Bool = false,
     )
 
-An AbstractNLPEvaluator for the problem:
-```
-min x1 * x4 * (x1 + x2 + x3) + x3
-st  x1 * x2 * x3 * x4 >= 25
-    x1^2 + x2^2 + x3^2 + x4^2 = 40
-    1 <= x1, x2, x3, x4 <= 5
+An [`MOI.AbstractNLPEvaluator`](@ref) for the problem:
+
+```math
+\\begin{aligned}
+\\text{min}       \\ & x_1 * x_4 * (x_1 + x_2 + x_3) + x_3 \\\\
+\\text{subject to}\\ & x_1 * x_2 * x_3 * x_4 \\ge 25 \\\\
+                     & x_1^2 + x_2^2 + x_3^2 + x_4^2 = 40 \\\\
+                     & 1 \\le x_1, x_2, x_3, x_4 \\le 5
+\\end{aligned}
 ```
-Start at (1,5,5,1)
-End at (1.000..., 4.743..., 3.821..., 1.379...)
+
+The optimal solution is `[1.000, 4.743, 3.821, 1.379]`.
 """
 struct HS071 <: MOI.AbstractNLPEvaluator
     enable_hessian::Bool
@@ -35,13 +38,13 @@ function MOI.initialize(d::HS071, requested_features::Vector{Symbol})
     for feat in requested_features
         if !(feat in MOI.features_available(d))
             error("Unsupported feature $feat")
-            # TODO: implement Jac-vec products for solvers that need them
         end
     end
+    return
 end
 
 function MOI.features_available(d::HS071)
-    features = [:Grad, :Jac, :ExprGraph]
+    features = [:Grad, :Jac, :JacVec, :ExprGraph]
     if d.enable_hessian
         push!(features, :Hess)
     end
@@ -99,6 +102,26 @@ function MOI.eval_objective_gradient(::HS071, grad_f, x)
     return
 end
 
+function MOI.constraint_gradient_structure(::HS071, ::Int)
+    return [1, 2, 3, 4]
+end
+
+function MOI.eval_constraint_gradient(::HS071, ∇g, x, i)
+    @assert 1 <= i <= 2
+    if i == 1
+        ∇g[1] = x[2] * x[3] * x[4]  # 1,1
+        ∇g[2] = x[1] * x[3] * x[4]  # 1,2
+        ∇g[3] = x[1] * x[2] * x[4]  # 1,3
+        ∇g[4] = x[1] * x[2] * x[3]  # 1,4
+    else
+        ∇g[1] = 2 * x[1]  # 2,1
+        ∇g[2] = 2 * x[2]  # 2,2
+        ∇g[3] = 2 * x[3]  # 2,3
+        ∇g[4] = 2 * x[4]  # 2,4
+    end
+    return
+end
+
 function MOI.jacobian_structure(::HS071)
     return Tuple{Int64,Int64}[
         (1, 1),
@@ -112,22 +135,6 @@ function MOI.jacobian_structure(::HS071)
     ]
 end
 
-function MOI.hessian_lagrangian_structure(d::HS071)
-    @assert d.enable_hessian
-    return Tuple{Int64,Int64}[
-        (1, 1),
-        (2, 1),
-        (2, 2),
-        (3, 1),
-        (3, 2),
-        (3, 3),
-        (4, 1),
-        (4, 2),
-        (4, 3),
-        (4, 4),
-    ]
-end
-
 function MOI.eval_constraint_jacobian(::HS071, J, x)
     # Constraint (row) 1
     J[1] = x[2] * x[3] * x[4]  # 1,1
@@ -142,6 +149,17 @@ function MOI.eval_constraint_jacobian(::HS071, J, x)
     return
 end
 
+function MOI.eval_constraint_jacobian_product(d::HS071, y, x, w)
+    y .= zero(eltype(y))
+    indices = MOI.jacobian_structure(d)
+    J = zeros(length(indices))
+    MOI.eval_constraint_jacobian(d, J, x)
+    for ((i, j), val) in zip(indices, J)
+        y[i] += val * w[j]
+    end
+    return
+end
+
 function MOI.eval_constraint_jacobian_transpose_product(::HS071, y, x, w)
     y[1] = (x[2] * x[3] * x[4]) * w[1] + (2 * x[1]) * w[2]
     y[2] = (x[1] * x[3] * x[4]) * w[1] + (2 * x[2]) * w[2]
@@ -150,6 +168,64 @@ function MOI.eval_constraint_jacobian_transpose_product(::HS071, y, x, w)
     return
 end
 
+function MOI.hessian_objective_structure(d::HS071)
+    return [(1, 1), (2, 1), (3, 1), (4, 1), (4, 2), (4, 3)]
+end
+
+function MOI.eval_hessian_objective(d::HS071, H, x)
+    @assert d.enable_hessian
+    H[1] = 2 * x[4]                 # 1,1
+    H[2] = x[4]                     # 2,1
+    H[3] = x[4]                     # 3,1
+    H[4] = 2 * x[1] + x[2] + x[3]   # 4,1
+    H[5] = x[1]                     # 4,2
+    H[6] = x[1]                     # 4,3
+    return
+end
+
+function MOI.hessian_constraint_structure(d::HS071, i::Int)
+    @assert 1 <= i <= 2
+    if i == 1
+        return [(2, 1), (3, 1), (3, 2), (4, 1), (4, 2), (4, 3)]
+    else
+        return [(1, 1), (2, 2), (3, 3), (4, 4)]
+    end
+end
+
+function MOI.eval_hessian_constraint(d::HS071, H, x, i)
+    @assert d.enable_hessian
+    if i == 1
+        H[1] = x[3] * x[4]  # 2,1
+        H[2] = x[2] * x[4]  # 3,1
+        H[3] = x[1] * x[4]  # 3,2
+        H[4] = x[2] * x[3]  # 4,1
+        H[5] = x[1] * x[3]  # 4,2
+        H[6] = x[1] * x[2]  # 4,3
+    else
+        H[1] = 2.0  # 1,1
+        H[2] = 2.0  # 2,2
+        H[3] = 2.0  # 3,3
+        H[4] = 2.0  # 4,4
+    end
+    return
+end
+
+function MOI.hessian_lagrangian_structure(d::HS071)
+    @assert d.enable_hessian
+    return Tuple{Int64,Int64}[
+        (1, 1),
+        (2, 1),
+        (2, 2),
+        (3, 1),
+        (3, 2),
+        (3, 3),
+        (4, 1),
+        (4, 2),
+        (4, 3),
+        (4, 4),
+    ]
+end
+
 function MOI.eval_hessian_lagrangian(d::HS071, H, x, σ, μ)
     @assert d.enable_hessian
     # Again, only lower left triangle
@@ -175,7 +251,7 @@ function MOI.eval_hessian_lagrangian(d::HS071, H, x, σ, μ)
     H[1] += μ[2] * 2  # 1,1
     H[3] += μ[2] * 2  # 2,2
     H[6] += μ[2] * 2  # 3,3
-    H[10] += μ[2] * 2
+    H[10] += μ[2] * 2 # 4,4
     return
 end
 
diff --git a/src/nlp.jl b/src/nlp.jl
index 255c3a18c3..31b328e577 100644
--- a/src/nlp.jl
+++ b/src/nlp.jl
@@ -19,6 +19,19 @@ Abstract supertype for the callback object that is used to query function
 values, derivatives, and expression graphs.
 
 It is used in [`NLPBlockData`](@ref).
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> supertype(typeof(evaluator))
+MathOptInterface.AbstractNLPEvaluator
+```
 """
 abstract type AbstractNLPEvaluator end
 
@@ -28,6 +41,24 @@ abstract type AbstractNLPEvaluator end
 An [`AbstractModelAttribute`](@ref) that stores an [`NLPBlockData`](@ref),
 representing a set of nonlinear constraints, and optionally a nonlinear
 objective.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}());
+
+julia> block = MOI.NLPBlockData(
+           MOI.NLPBoundsPair.([25.0, 40.0], [Inf, 40.0]),
+           MOI.Test.HS071(true),
+           true,
+       );
+
+julia> MOI.set(model, MOI.NLPBlock(), block)
+```
 """
 struct NLPBlock <: AbstractModelAttribute end
 
@@ -38,6 +69,18 @@ An [`AbstractModelAttribute`](@ref) for the Lagrange multipliers on the
 constraints from the [`NLPBlock`](@ref) in result `result_index`.
 
 If `result_index` is omitted, it is `1` by default.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> MOI.NLPBlockDual()
+MathOptInterface.NLPBlockDual(1)
+
+julia> MOI.NLPBlockDual(2)
+MathOptInterface.NLPBlockDual(2)
+```
 """
 struct NLPBlockDual <: AbstractModelAttribute
     result_index::Int
@@ -52,6 +95,26 @@ is_set_by_optimize(::NLPBlockDual) = true
 An [`AbstractModelAttribute`](@ref) for the initial assignment of the Lagrange
 multipliers on the constraints from the [`NLPBlock`](@ref) that the solver may
 use to warm-start the solve.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}());
+
+julia> block = MOI.NLPBlockData(
+           MOI.NLPBoundsPair.([25.0, 40.0], [Inf, 40.0]),
+           MOI.Test.HS071(true),
+           true,
+       );
+
+julia> MOI.set(model, MOI.NLPBlock(), block)
+
+julia> MOI.set(model, MOI.NLPBlockDualStart(), [1.0, 2.0])
+```
 """
 struct NLPBlockDualStart <: AbstractModelAttribute end
 
@@ -61,6 +124,17 @@ struct NLPBlockDualStart <: AbstractModelAttribute end
 A struct holding a pair of lower and upper bounds.
 
 `-Inf` and `Inf` can be used to indicate no lower or upper bound, respectively.
+
+## Example
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> bounds = MOI.NLPBoundsPair.([25.0, 40.0], [Inf, 40.0])
+2-element Vector{MathOptInterface.NLPBoundsPair}:
+ MathOptInterface.NLPBoundsPair(25.0, Inf)
+ MathOptInterface.NLPBoundsPair(40.0, 40.0)
+```
 """
 struct NLPBoundsPair
     lower::Float64
@@ -90,6 +164,24 @@ Hessian-of-the-Lagrangian queries, `σ` must be set to zero.
     Throughout the evaluator, all variables are ordered according to
     [`ListOfVariableIndices`](@ref). Hence, MOI copies of nonlinear problems
     must not re-order variables.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> model = MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}());
+
+julia> block = MOI.NLPBlockData(
+           MOI.NLPBoundsPair.([25.0, 40.0], [Inf, 40.0]),
+           MOI.Test.HS071(true),
+           true,
+       );
+
+julia> MOI.set(model, MOI.NLPBlock(), block)
+```
 """
 struct NLPBlockData
     constraint_bounds::Vector{NLPBoundsPair}
@@ -115,7 +207,8 @@ are supported by `d`.
 The following features are defined:
 
  * `:Grad`: enables [`eval_objective_gradient`](@ref)
- * `:Jac`: enables [`eval_constraint_jacobian`](@ref)
+ * `:Jac`: enables [`eval_constraint_jacobian`](@ref) and
+   [`eval_constraint_gradient`](@ref)
  * `:JacVec`: enables [`eval_constraint_jacobian_product`](@ref) and
    [`eval_constraint_jacobian_transpose_product`](@ref)
  * `:Hess`: enables [`eval_hessian_lagrangian`](@ref)
@@ -125,12 +218,16 @@ The following features are defined:
 In all cases, including when `requested_features` is empty,
 [`eval_objective`](@ref) and [`eval_constraint`](@ref) are supported.
 
-## Examples
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
 
-```julia
-MOI.initialize(d, Symbol[])
-MOI.initialize(d, [:ExprGraph])
-MOI.initialize(d, MOI.features_available(d))
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, [:Grad, :Jac])
 ```
 """
 function initialize end
@@ -141,6 +238,25 @@ function initialize end
 Returns the subset of features available for this problem instance.
 
 See [`initialize`](@ref) for the list of defined features.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true, true);
+
+julia> MOI.features_available(evaluator)
+6-element Vector{Symbol}:
+ :Grad
+ :Jac
+ :JacVec
+ :ExprGraph
+ :Hess
+ :HessVec
+```
 """
 function features_available end
 
@@ -148,11 +264,32 @@ function features_available end
     eval_objective(d::AbstractNLPEvaluator, x::AbstractVector{T})::T where {T}
 
 Evaluate the objective ``f(x)``, returning a scalar value.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref), but you do not
+need to pass a value.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[])
+
+julia> MOI.eval_objective(evaluator, [1.0, 2.0, 3.0, 4.0])
+27.0
+```
 """
 function eval_objective end
 
 """
-    eval_constraint(d::AbstractNLPEvaluator,
+    eval_constraint(
+        d::AbstractNLPEvaluator,
         g::AbstractVector{T},
         x::AbstractVector{T},
     )::Nothing where {T}
@@ -160,11 +297,37 @@ function eval_objective end
 Given a set of vector-valued constraints ``l \\le g(x) \\le u``, evaluate the
 constraint function ``g(x)``, storing the result in the vector `g`.
 
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref), but you do not
+need to pass a value.
+
 ## Implementation notes
 
 When implementing this method, you must not assume that `g` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[])
+
+julia> g = fill(NaN, 2);
+
+julia> MOI.eval_constraint(evaluator, g, [1.0, 2.0, 3.0, 4.0])
+
+julia> g
+2-element Vector{Float64}:
+ 24.0
+ 30.0
+```
 """
 function eval_constraint end
 
@@ -178,11 +341,38 @@ function eval_constraint end
 Evaluate the gradient of the objective function ``grad = \\nabla f(x)`` as a
 dense vector, storing the result in the vector `grad`.
 
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Grad`.
+
 ## Implementation notes
 
 When implementing this method, you must not assume that `grad` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Grad])
+
+julia> grad = fill(NaN, 4);
+
+julia> MOI.eval_objective_gradient(evaluator, grad, [1.0, 2.0, 3.0, 4.0])
+
+julia> grad
+4-element Vector{Float64}:
+ 28.0
+  4.0
+  5.0
+  6.0
+```
 """
 function eval_objective_gradient end
 
@@ -200,6 +390,33 @@ case the solver should combine the corresponding elements by adding them
 together.
 
 The sparsity structure is assumed to be independent of the point ``x``.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Jac`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac])
+
+julia> MOI.jacobian_structure(evaluator)
+8-element Vector{Tuple{Int64, Int64}}:
+ (1, 1)
+ (1, 2)
+ (1, 3)
+ (1, 4)
+ (2, 1)
+ (2, 2)
+ (2, 3)
+ (2, 4)
+```
 """
 function jacobian_structure end
 
@@ -216,10 +433,35 @@ The indices are not required to be sorted and can contain duplicates, in which
 case the solver should combine the corresponding elements by adding them
 together.
 
-Any mix of lower and upper-triangular indices is valid. Elements `(i,j)` and
-`(j,i)`, if both present, should be treated as duplicates.
+Any mix of lower and upper-triangular indices is valid. Elements `(i, j)` and
+`(j, i)`, if both present, should be treated as duplicates.
 
 The sparsity structure is assumed to be independent of the point ``x``.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> MOI.hessian_objective_structure(evaluator)
+6-element Vector{Tuple{Int64, Int64}}:
+ (1, 1)
+ (2, 1)
+ (3, 1)
+ (4, 1)
+ (4, 2)
+ (4, 3)
+```
 """
 function hessian_objective_structure end
 
@@ -237,10 +479,42 @@ The indices are not required to be sorted and can contain duplicates, in which
 case the solver should combine the corresponding elements by adding them
 together.
 
-Any mix of lower and upper-triangular indices is valid. Elements `(i,j)` and
-`(j,i)`, if both present, should be treated as duplicates.
+Any mix of lower and upper-triangular indices is valid. Elements `(i, j)` and
+`(j, i)`, if both present, should be treated as duplicates.
 
 The sparsity structure is assumed to be independent of the point ``x``.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> MOI.hessian_constraint_structure(evaluator, 1)
+6-element Vector{Tuple{Int64, Int64}}:
+ (2, 1)
+ (3, 1)
+ (3, 2)
+ (4, 1)
+ (4, 2)
+ (4, 3)
+
+julia> MOI.hessian_constraint_structure(evaluator, 2)
+4-element Vector{Tuple{Int64, Int64}}:
+ (1, 1)
+ (2, 2)
+ (3, 3)
+ (4, 4)
+```
 """
 function hessian_constraint_structure end
 
@@ -257,10 +531,39 @@ The indices are not required to be sorted and can contain duplicates, in which
 case the solver should combine the corresponding elements by adding them
 together.
 
-Any mix of lower and upper-triangular indices is valid. Elements `(i,j)` and
-`(j,i)`, if both present, should be treated as duplicates.
+Any mix of lower and upper-triangular indices is valid. Elements `(i, j)` and
+`(j, i)`, if both present, should be treated as duplicates.
 
 The sparsity structure is assumed to be independent of the point ``x``.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> MOI.hessian_lagrangian_structure(evaluator)
+10-element Vector{Tuple{Int64, Int64}}:
+ (1, 1)
+ (2, 1)
+ (2, 2)
+ (3, 1)
+ (3, 2)
+ (3, 3)
+ (4, 1)
+ (4, 2)
+ (4, 3)
+ (4, 4)
+```
 """
 function hessian_lagrangian_structure end
 
@@ -275,6 +578,29 @@ case the solver should combine the corresponding elements by adding them
 together.
 
 The sparsity structure is assumed to be independent of the point ``x``.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Jac`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac])
+
+julia> indices = MOI.constraint_gradient_structure(evaluator, 1)
+4-element Vector{Int64}:
+ 1
+ 2
+ 3
+ 4
+```
 """
 function constraint_gradient_structure end
 
@@ -295,11 +621,41 @@ non-zero values in `∇g`, corresponding to the structure returned by
 When implementing this method, you must not assume that `∇g` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Jac`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac])
+
+julia> indices = MOI.constraint_gradient_structure(evaluator, 1);
+
+julia> ∇g = zeros(length(indices));
+
+julia> MOI.eval_constraint_gradient(evaluator, ∇g, [1.0, 2.0, 3.0, 4.0], 1)
+
+julia> ∇g
+4-element Vector{Float64}:
+ 24.0
+ 12.0
+  8.0
+  6.0
+```
 """
 function eval_constraint_gradient end
 
 """
-    eval_constraint_jacobian(d::AbstractNLPEvaluator,
+    eval_constraint_jacobian(
+        d::AbstractNLPEvaluator,
         J::AbstractVector{T},
         x::AbstractVector{T},
     )::Nothing where {T}
@@ -315,6 +671,39 @@ by [`jacobian_structure`](@ref).
 When implementing this method, you must not assume that `J` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac])
+
+julia> J_indices = MOI.jacobian_structure(evaluator);
+
+julia> J = zeros(length(J_indices));
+
+julia> MOI.eval_constraint_jacobian(evaluator, J, [1.0, 2.0, 3.0, 4.0])
+
+julia> J
+8-element Vector{Float64}:
+ 24.0
+ 12.0
+  8.0
+  6.0
+  2.0
+  4.0
+  6.0
+  8.0
+```
 """
 function eval_constraint_jacobian end
 
@@ -337,6 +726,35 @@ is the number of nonlinear constraints.
 When implementing this method, you must not assume that `y` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:JacVec`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac, :JacVec])
+
+julia> y = zeros(2);
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> w = [1.5, 2.5, 3.5, 4.5];
+
+julia> MOI.eval_constraint_jacobian_product(evaluator, y, x, w)
+
+julia> y
+2-element Vector{Float64}:
+ 121.0
+  70.0
+```
 """
 function eval_constraint_jacobian_product end
 
@@ -359,6 +777,37 @@ is the number of nonlinear constraints.
 When implementing this method, you must not assume that `y` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:JacVec`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Jac, :JacVec])
+
+julia> y = zeros(4);
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> w = [1.5, 2.5];
+
+julia> MOI.eval_constraint_jacobian_transpose_product(evaluator, y, x, w)
+
+julia> y
+4-element Vector{Float64}:
+ 41.0
+ 28.0
+ 27.0
+ 29.0
+```
 """
 function eval_constraint_jacobian_transpose_product end
 
@@ -384,6 +833,41 @@ The vectors have dimensions such that `length(h) == length(x) == length(v)`.
 When implementing this method, you must not assume that `h` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:HessVec`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true, true);
+
+julia> MOI.initialize(evaluator, Symbol[:HessVec])
+
+julia> H = fill(NaN, 4);
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> v = [1.5, 2.5, 3.5, 4.5];
+
+julia> σ = 1.0;
+
+julia> μ = [1.0, 1.0];
+
+julia> MOI.eval_hessian_lagrangian_product(evaluator, H, x, v, σ, μ)
+
+julia> H
+4-element Vector{Float64}:
+ 155.5
+  61.0
+  48.5
+  49.0
+```
 """
 function eval_hessian_lagrangian_product end
 
@@ -404,6 +888,39 @@ returned by [`hessian_objective_structure`](@ref).
 When implementing this method, you must not assume that `H` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true, true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> indices = MOI.hessian_objective_structure(evaluator);
+
+julia> H = zeros(length(indices));
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> MOI.eval_hessian_objective(evaluator, H, x)
+
+julia> H
+6-element Vector{Float64}:
+ 8.0
+ 4.0
+ 4.0
+ 7.0
+ 1.0
+ 1.0
+```
 """
 function eval_hessian_objective end
 
@@ -425,6 +942,39 @@ returned by [`hessian_constraint_structure`](@ref).
 When implementing this method, you must not assume that `H` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true, true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> indices = MOI.hessian_constraint_structure(evaluator, 1);
+
+julia> H = zeros(length(indices));
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> MOI.eval_hessian_constraint(evaluator, H, x, 1)
+
+julia> H
+6-element Vector{Float64}:
+ 12.0
+  8.0
+  4.0
+  6.0
+  3.0
+  2.0
+```
 """
 function eval_hessian_constraint end
 
@@ -448,6 +998,47 @@ returned by [`hessian_lagrangian_structure`](@ref).
 When implementing this method, you must not assume that `H` is
 `Vector{Float64}`, but you may assume that it supports `setindex!` and `length`.
 For example, it may be the `view` of a vector.
+
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with `:Hess`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, Symbol[:Hess])
+
+julia> indices = MOI.hessian_lagrangian_structure(evaluator);
+
+julia> H = zeros(length(indices));
+
+julia> x = [1.0, 2.0, 3.0, 4.0];
+
+julia> σ = 1.0;
+
+julia> μ = [1.0, 1.0];
+
+julia> MOI.eval_hessian_lagrangian(evaluator, H, x, σ, μ)
+
+julia> H
+10-element Vector{Float64}:
+ 10.0
+ 16.0
+  2.0
+ 12.0
+  4.0
+  2.0
+ 13.0
+  4.0
+  3.0
+  2.0
+```
 """
 function eval_hessian_lagrangian end
 
@@ -472,25 +1063,24 @@ expressions:
    etc., but modeling interfaces may choose to extend these basic functions, or
    error if they encounter unsupported functions.
 
-## Examples
+## Initialize
 
-The expression ``x_1+\\sin(x_2/\\exp(x_3))`` is represented as
-```julia
-:(x[MOI.VariableIndex(1)] + sin(x[MOI.VariableIndex(2)] / exp(x[MOI.VariableIndex[3]])))
-```
-or equivalently
-```julia
-Expr(
-    :call,
-    :+,
-    Expr(:ref, :x, MOI.VariableIndex(1)),
-    Expr(
-        :call,
-        :/,
-        Expr(:call, :sin, Expr(:ref, :x, MOI.VariableIndex(2))),
-        Expr(:call, :exp, Expr(:ref, :x, MOI.VariableIndex(3))),
-    ),
-)
+Before calling this function, you must call [`initialize`](@ref) with
+`:ExprGraph`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, [:ExprGraph])
+
+julia> MOI.objective_expr(evaluator)
+:(x[MOI.VariableIndex(1)] * x[MOI.VariableIndex(4)] * (x[MOI.VariableIndex(1)] + x[MOI.VariableIndex(2)] + x[MOI.VariableIndex(3)]) + x[MOI.VariableIndex(3)])
 ```
 """
 function objective_expr end
@@ -509,20 +1099,34 @@ comparison operator indicating the sense of and bounds on the constraint.
 For single-sided comparisons, the body of the constraint must be on the
 left-hand side, and the right-hand side must be a constant.
 
-For double-sided comparisons (that is, ``l \\le f(x) \\le u``), the body of the
+For double-sided comparisons (that is, ``l \\le g(x) \\le u``), the body of the
 constraint must be in the middle, and the left- and right-hand sides must be
 constants.
 
 The bounds on the constraints must match the [`NLPBoundsPair`](@ref)s passed to
 [`NLPBlockData`](@ref).
 
-## Examples
+## Initialize
+
+Before calling this function, you must call [`initialize`](@ref) with
+`:ExprGraph`.
+
+## Example
+
+This example uses the [`Test.HS071`](@ref) evaluator.
+
+```jldoctest
+julia> import MathOptInterface as MOI
+
+julia> evaluator = MOI.Test.HS071(true);
+
+julia> MOI.initialize(evaluator, [:ExprGraph])
+
+julia> MOI.constraint_expr(evaluator, 1)
+:(x[MOI.VariableIndex(1)] * x[MOI.VariableIndex(2)] * x[MOI.VariableIndex(3)] * x[MOI.VariableIndex(4)] >= 25.0)
 
-```julia
-:(x[MOI.VariableIndex(1)]^2 <= 1.0)
-:(x[MOI.VariableIndex(1)]^2 >= 2.0)
-:(x[MOI.VariableIndex(1)]^2 == 3.0)
-:(4.0 <= x[MOI.VariableIndex(1)]^2 <= 5.0)
+julia> MOI.constraint_expr(evaluator, 2)
+:(x[MOI.VariableIndex(1)] ^ 2 + x[MOI.VariableIndex(2)] ^ 2 + x[MOI.VariableIndex(3)] ^ 2 + x[MOI.VariableIndex(4)] ^ 2 == 40.0)
 ```
 """
 function constraint_expr end
diff --git a/test/Test/Test.jl b/test/Test/Test.jl
index fbadb3a6a4..d3238df04a 100644
--- a/test/Test/Test.jl
+++ b/test/Test/Test.jl
@@ -124,3 +124,75 @@ end
     exclude = [r"^test_model_Name$", r"test_.+"]
     @test_logs MOI.Test.runtests(model, config; exclude = exclude)
 end
+
+@testset "test_HS071_evaluator" begin
+    evaluator = MOI.Test.HS071(true, true)
+    features = [:Grad, :Jac, :JacVec, :ExprGraph, :Hess, :HessVec]
+    @test MOI.features_available(evaluator) == features
+    @test_throws(
+        ErrorException("Unsupported feature foo"),
+        MOI.initialize(evaluator, [:foo]),
+    )
+    MOI.initialize(evaluator, features)
+    x = [1.0, 2.0, 3.0, 4.0]
+    @test MOI.eval_objective(evaluator, x) == 27.0
+    g = fill(NaN, 2)
+    MOI.eval_constraint(evaluator, g, x)
+    @test g == [24.0, 30.0]
+    grad = fill(NaN, 4)
+    MOI.eval_objective_gradient(evaluator, grad, x)
+    @test grad == [28.0, 4.0, 5.0, 6.0]
+    @test MOI.jacobian_structure(evaluator) ==
+          [(1, 1), (1, 2), (1, 3), (1, 4), (2, 1), (2, 2), (2, 3), (2, 4)]
+    @test MOI.hessian_objective_structure(evaluator) ==
+          [(1, 1), (2, 1), (3, 1), (4, 1), (4, 2), (4, 3)]
+    @test MOI.hessian_constraint_structure(evaluator, 1) ==
+          [(2, 1), (3, 1), (3, 2), (4, 1), (4, 2), (4, 3)]
+    @test MOI.hessian_constraint_structure(evaluator, 2) ==
+          [(1, 1), (2, 2), (3, 3), (4, 4)]
+    @test MOI.hessian_lagrangian_structure(evaluator) == vcat(
+        [(1, 1), (2, 1), (2, 2), (3, 1), (3, 2), (3, 3)],
+        [(4, 1), (4, 2), (4, 3), (4, 4)],
+    )
+    @test MOI.constraint_gradient_structure(evaluator, 1) == [1, 2, 3, 4]
+    @test MOI.constraint_gradient_structure(evaluator, 2) == [1, 2, 3, 4]
+    ∇g = zeros(4)
+    MOI.eval_constraint_gradient(evaluator, ∇g, x, 1)
+    @test ∇g == [24.0, 12.0, 8.0, 6.0]
+    MOI.eval_constraint_gradient(evaluator, ∇g, x, 2)
+    @test ∇g == [2.0, 4.0, 6.0, 8.0]
+    J = zeros(8)
+    MOI.eval_constraint_jacobian(evaluator, J, x)
+    @test J == [24.0, 12.0, 8.0, 6.0, 2.0, 4.0, 6.0, 8.0]
+    y = zeros(2)
+    w = [1.5, 2.5, 3.5, 4.5]
+    MOI.eval_constraint_jacobian_product(evaluator, y, x, w)
+    @test y == [121.0, 70.0]
+    y = zeros(4)
+    w2 = @view(w[1:2])
+    MOI.eval_constraint_jacobian_transpose_product(evaluator, y, x, w2)
+    @test y == [41.0, 28.0, 27.0, 29.0]
+    σ = 1.0
+    μ = [1.0, 1.0]
+    H = zeros(4)
+    MOI.eval_hessian_lagrangian_product(evaluator, H, x, w, σ, μ)
+    @test H == [155.5, 61.0, 48.5, 49.0]
+    H = zeros(6)
+    MOI.eval_hessian_objective(evaluator, H, x)
+    @test H == [8.0, 4.0, 4.0, 7.0, 1.0, 1.0]
+    MOI.eval_hessian_constraint(evaluator, H, x, 1)
+    @test H == [12.0, 8.0, 4.0, 6.0, 3.0, 2.0]
+    H = zeros(4)
+    MOI.eval_hessian_constraint(evaluator, H, x, 2)
+    @test H == [2.0, 2.0, 2.0, 2.0]
+    H = zeros(10)
+    MOI.eval_hessian_lagrangian(evaluator, H, x, σ, μ)
+    @test H == [10.0, 16.0, 2.0, 12.0, 4.0, 2.0, 13.0, 4.0, 3.0, 2.0]
+    x = Expr.(:ref, :x, MOI.VariableIndex.(1:4))
+    @test MOI.objective_expr(evaluator) ==
+          :($(x[1]) * $(x[4]) * ($(x[1]) + $(x[2]) + $(x[3])) + $(x[3]))
+    @test MOI.constraint_expr(evaluator, 1) ==
+          :($(x[1]) * $(x[2]) * $(x[3]) * $(x[4]) >= 25.0)
+    @test MOI.constraint_expr(evaluator, 2) ==
+          :($(x[1])^2 + $(x[2])^2 + $(x[3])^2 + $(x[4])^2 == 40.0)
+end