diff --git a/README.md b/README.md
index 3557944b4..a6781c55c 100644
--- a/README.md
+++ b/README.md
@@ -70,7 +70,7 @@ If Jacobian-vector products can be computed more efficiently than by evaluating
 
 The following method is defined if second-order derivatives are available:
 
-* `hess(model, x, y)`: evaluate *∇²L(x,y)*, the Hessian of the Lagrangian at `x` and `y`
+* `hess(model, x, y)`: evaluate *∇²L(x,y)*, the Hessian of the Lagrangian at `x` and `y` as a sparse matrix
 
 If Hessian-vector products can be computed more efficiently than by evaluating the Hessian explicitly, the following method may be implemented:
 
@@ -111,16 +111,18 @@ Attribute          | Type          | Notes
 `jfree`            | `Vector{Int}` | indices of "free" constraints (there shouldn't be any)
 `jinf`             | `Vector{Int}` | indices of the visibly infeasible constraints
 `nnzo`             | `Int`         | number of nonzeros in the gradient
-`nnzj`             | `Int`         | number of nonzeros in the sparse Jacobian
-`lin_nnzj`         | `Int`         | number of nonzeros in the sparse linear constraints Jacobian
-`nln_nnzj`         | `Int`         | number of nonzeros in the sparse nonlinear constraints Jacobian
-`nnzh`             | `Int`         | number of nonzeros in the lower triangular part of the sparse Hessian of the Lagrangian
+`nnzj`             | `Int`         | number of nonzeros in the Jacobian
+`lin_nnzj`         | `Int`         | number of nonzeros in the linear constraints Jacobian
+`nln_nnzj`         | `Int`         | number of nonzeros in the nonlinear constraints Jacobian
+`nnzh`             | `Int`         | number of nonzeros in the lower triangular part of the Hessian of the Lagrangian
 `minimize`         | `Bool`        | true if `optimize == minimize`
 `islp`             | `Bool`        | true if the problem is a linear program
 `name`             | `String`      | problem name
+`sparse_jacobian`  | `Bool`        | true if the Jacobian of the constraints is sparse
+`sparse_hessian`   | `Bool`        | true if the Hessian of the Lagrangian is sparse
 `grad_available`   | `Bool`        | true if the gradient of the objective is available
-`jac_available`    | `Bool`        | true if the sparse Jacobian of the constraints is available
-`hess_available`   | `Bool`        | true if the sparse Hessian of the Lagrangian is available
+`jac_available`    | `Bool`        | true if the Jacobian of the constraints is available
+`hess_available`   | `Bool`        | true if the Hessian of the Lagrangian is available
 `jprod_available`  | `Bool`        | true if the Jacobian-vector product `J * v` is available
 `jtprod_available` | `Bool`        | true if the transpose Jacobian-vector product `J' * v` is available
 `hprod_available`  | `Bool`        | true if the Hessian-vector product of the Lagrangian `H * v` is available
diff --git a/docs/src/api.md b/docs/src/api.md
index c7424fb87..cd91de806 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -29,6 +29,7 @@ If not, click on the link and go to the description.
 
 - `!` means inplace;
 - `_coord` means coordinate format;
+- `_dense` means dense format;
 - `prod` means matrix-vector product;
 - `_op` means operator (as in [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl));
 - `_lin` and `_nln` respectively refer to linear and nonlinear constraints.
@@ -40,13 +41,14 @@ NLPModels instances.
 |-------------------|-------------------------------------------|
 | ``f(x)``            | [`obj`](@ref), [`objgrad`](@ref), [`objgrad!`](@ref), [`objcons`](@ref), [`objcons!`](@ref) |
 | ``\nabla f(x)``     | [`grad`](@ref), [`grad!`](@ref), [`objgrad`](@ref), [`objgrad!`](@ref) |
-| ``\nabla^2 f(x)``   | [`hess`](@ref), [`hess_op`](@ref), [`hess_op!`](@ref), [`hess_coord`](@ref), [`hess_coord`](@ref), [`hess_structure`](@ref), [`hess_structure!`](@ref), [`hprod`](@ref), [`hprod!`](@ref) |
+| ``\nabla^2 f(x)``   | [`hess`](@ref), [`hess_op`](@ref), [`hess_op!`](@ref), [`hess_coord`](@ref), [`hess_coord!`](@ref), [`hess_dense!`](@ref), [`hess_structure`](@ref), [`hess_structure!`](@ref), [`hprod`](@ref), [`hprod!`](@ref) |
 | ``c(x)``            | [`cons_lin`](@ref), [`cons_lin!`](@ref), [`cons_nln`](@ref), [`cons_nln!`](@ref), [`cons`](@ref), [`cons!`](@ref), [`objcons`](@ref), [`objcons!`](@ref) |
-| ``J(x)``            | [`jac_lin`](@ref), [`jac_nln`](@ref), [`jac`](@ref), [`jac_lin_op`](@ref), [`jac_lin_op!`](@ref), [`jac_nln_op`](@ref), [`jac_nln_op!`](@ref),[`jac_op`](@ref), [`jac_op!`](@ref), [`jac_lin_coord`](@ref), [`jac_lin_coord!`](@ref), [`jac_nln_coord`](@ref), [`jac_nln_coord!`](@ref), [`jac_coord`](@ref), [`jac_coord!`](@ref), [`jac_lin_structure`](@ref), [`jac_lin_structure!`](@ref), [`jac_nln_structure`](@ref), [`jac_nln_structure!`](@ref), [`jac_structure`](@ref), [`jprod_lin`](@ref), [`jprod_lin!`](@ref), [`jprod_nln`](@ref), [`jprod_nln!`](@ref), [`jprod`](@ref), [`jprod!`](@ref), [`jtprod_lin`](@ref), [`jtprod_lin!`](@ref), [`jtprod_nln`](@ref), [`jtprod_nln!`](@ref), [`jtprod`](@ref), [`jtprod!`](@ref) |
-| ``\nabla^2 L(x,y)`` | [`hess`](@ref), [`hess_op`](@ref), [`hess_coord`](@ref), [`hess_coord!`](@ref), [`hess_structure`](@ref), [`hess_structure!`](@ref), [`hprod`](@ref), [`hprod!`](@ref), [`jth_hprod`](@ref), [`jth_hprod!`](@ref), [`jth_hess`](@ref), [`jth_hess_coord`](@ref), [`jth_hess_coord!`](@ref), [`ghjvprod`](@ref), [`ghjvprod!`](@ref)  |
+| ``J(x)``            | [`jac_lin`](@ref), [`jac_nln`](@ref), [`jac`](@ref), [`jac_lin_op`](@ref), [`jac_lin_op!`](@ref), [`jac_nln_op`](@ref), [`jac_nln_op!`](@ref),[`jac_op`](@ref), [`jac_op!`](@ref), [`jac_lin_coord`](@ref), [`jac_lin_coord!`](@ref), [`jac_nln_coord`](@ref), [`jac_nln_coord!`](@ref), [`jac_coord`](@ref), [`jac_coord!`](@ref), [`jac_dense!`](@ref), [`jac_lin_structure`](@ref), [`jac_lin_structure!`](@ref), [`jac_nln_structure`](@ref), [`jac_nln_structure!`](@ref), [`jac_structure`](@ref), [`jprod_lin`](@ref), [`jprod_lin!`](@ref), [`jprod_nln`](@ref), [`jprod_nln!`](@ref), [`jprod`](@ref), [`jprod!`](@ref), [`jtprod_lin`](@ref), [`jtprod_lin!`](@ref), [`jtprod_nln`](@ref), [`jtprod_nln!`](@ref), [`jtprod`](@ref), [`jtprod!`](@ref) |
+| ``\nabla^2 L(x,y)`` | [`hess`](@ref), [`hess_op`](@ref), [`hess_coord`](@ref), [`hess_coord!`](@ref), [`hess_dense!`](@ref), [`hess_structure`](@ref), [`hess_structure!`](@ref), [`hprod`](@ref), [`hprod!`](@ref), [`jth_hprod`](@ref), [`jth_hprod!`](@ref), [`jth_hess`](@ref), [`jth_hess_coord`](@ref), [`jth_hess_coord!`](@ref), [`ghjvprod`](@ref), [`ghjvprod!`](@ref)  |
 
 If only a subset of the functions listed above is implemented, you can indicate which ones are not available when creating the [`NLPModelMeta`](@ref), using the keyword arguments
 `grad_available`, `jac_available`, `hess_available`, `jprod_available`, `jtprod_available`, and `hprod_available`.
+You can also specify whether the Jacobian of the constraints and the Hessian of the objective or Lagrangian are sparse using the keyword arguments `sparse_jacobian` and `sparse_hessian`.
 
 ## [API for NLSModels](@id nls-api)
 
diff --git a/docs/src/guidelines.md b/docs/src/guidelines.md
index b5e426345..92b9c8d05 100644
--- a/docs/src/guidelines.md
+++ b/docs/src/guidelines.md
@@ -60,27 +60,42 @@ The following functions should be defined:
 - Objective (unconstrained models only need to worry about these)
   - `obj(nlp, x)`
   - `grad!(nlp, x, g)`
-  - `hess_structure!(nlp, hrows, hcols)`
-  - `hess_coord!(nlp, x, hvals; obj_weight=1)`
+  - `hess_structure!(nlp, hrows, hcols)` (sparse Hessian)
+  - `hess_coord!(nlp, x, hvals; obj_weight=1)` (sparse Hessian)
+  - `hess_dense!(nlp, x, Hx; obj_weight=1)` (dense Hessian)
   - `hprod!(nlp, x, v, Hv; obj_weight=1)` (actually defaults to calling the constrained case)
+
 - Constraints (constrained models need to worry about these and the ones above)
-  - `cons_lin!(nlp, x, c)`
-  - `cons_nln!(nlp, x, c)`
-  - `jac_lin_structure!(nlp, jrows, jcols)`
-  - `jac_nln_structure!(nlp, jrows, jcols)`
-  - `jac_lin_coord!(nlp, x, jvals)`
-  - `jac_nln_coord!(nlp, x, jvals)`
-  - `jprod_lin!(nlp, x, v, Jv)`
-  - `jprod_nln!(nlp, x, v, Jv)`
-  - `jtprod_lin!(nlp, x, v, Jtv)`
-  - `jtprod_nln!(nlp, x, v, Jtv)`
-  - `hess_coord!(nlp, x, y, hvals; obj_weight=1)`
+  - `cons!(nlp, x, c)`
+  - `jac_structure!(nlp, jrows, jcols)` (sparse Jacobian)
+  - `jac_coord!(nlp, x, jvals)` (sparse Jacobian)
+  - `jac_dense!(nlp, x, Jx)` (dense Jacobian)
+  - `jprod!(nlp, x, v, Jv)`
+  - `jtprod!(nlp, x, v, Jtv)`
+  - `hess_coord!(nlp, x, y, hvals; obj_weight=1)` (sparse Hessian)
+  - `hess_dense!(nlp, x, y, Hx; obj_weight=1)` (dense Hessian)
   - `hprod!(nlp, x, y, v, Hv; obj_weight=1)`
 
+When the split API is enabled (`nlp.meta.split_api = true`), the treatment of linear and nonlinear constraints and their associated Jacobians is separated, and the following additional functions must be provided:
+
+- `cons_lin!(nlp, x, c)`
+- `cons_nln!(nlp, x, c)`
+- `jac_lin_structure!(nlp, jrows, jcols)`
+- `jac_nln_structure!(nlp, jrows, jcols)`
+- `jac_lin_coord!(nlp, x, jvals)`
+- `jac_nln_coord!(nlp, x, jvals)`
+- `jprod_lin!(nlp, x, v, Jv)`
+- `jprod_nln!(nlp, x, v, Jv)`
+- `jtprod_lin!(nlp, x, v, Jtv)`
+- `jtprod_nln!(nlp, x, v, Jtv)`
+
 The linear constraints are specified at the initialization of the `NLPModelMeta` using the keyword arguement `lin`.
 The indices of linear and nonlinear constraints are respectively available in `nlp.meta.lin` and `nlp.meta.nln`.
-If your model uses only linear (resp. nonlinear) constraints, then it suffices to implement the `*_lin` (resp. `*_nln`) functions.
-Alternatively, one could implement only the functions without the suffixes `_nln!` (e.g., only `cons!`), but this might run into errors with tools differentiating linear and nonlinear constraints.
+
+If the Jacobian or the Hessian of the Lagrangian is dense, there is no need to implement the corresponding `*_structure!` and `*_coord!` methods.
+Only the corresponding `*_dense!` methods need to be implemented.
+Note that only the in-place API is provided for dense Jacobians and Hessians.
+This is specified at the initialization of [`NLPModelMeta`](@ref) through the keyword arguments `sparse_jacobian` and `sparse_hessian`.
 
 ## [Availability of the API](@id availability-api)
 
diff --git a/docs/src/index.md b/docs/src/index.md
index c5c4f123a..3a8758cc2 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -106,16 +106,18 @@ Attribute          | Type          | Notes
 `jfree`            | `Vector{Int}` | indices of "free" constraints (there shouldn't be any)
 `jinf`             | `Vector{Int}` | indices of the visibly infeasible constraints
 `nnzo`             | `Int`         | number of nonzeros in the gradient
-`nnzj`             | `Int`         | number of nonzeros in the sparse Jacobian
-`lin_nnzj`         | `Int`         | number of nonzeros in the sparse linear constraints Jacobian
-`nln_nnzj`         | `Int`         | number of nonzeros in the sparse nonlinear constraints Jacobian
-`nnzh`             | `Int`         | number of nonzeros in the lower triangular part of the sparse Hessian of the Lagrangian
+`nnzj`             | `Int`         | number of nonzeros in the Jacobian
+`lin_nnzj`         | `Int`         | number of nonzeros in the linear constraints Jacobian
+`nln_nnzj`         | `Int`         | number of nonzeros in the nonlinear constraints Jacobian
+`nnzh`             | `Int`         | number of nonzeros in the lower triangular part of the Hessian of the Lagrangian
 `minimize`         | `Bool`        | true if `optimize == minimize`
 `islp`             | `Bool`        | true if the problem is a linear program
 `name`             | `String`      | problem name
+`sparse_jacobian`  | `Bool`        | true if the Jacobian of the constraints is sparse
+`sparse_hessian`   | `Bool`        | true if the Hessian of the Lagrangian is sparse
 `grad_available`   | `Bool`        | true if the gradient of the objective is available
-`jac_available`    | `Bool`        | true if the sparse Jacobian of the constraints is available
-`hess_available`   | `Bool`        | true if the sparse Hessian of the Lagrangian is available
+`jac_available`    | `Bool`        | true if the Jacobian of the constraints is available
+`hess_available`   | `Bool`        | true if the Hessian of the Lagrangian is available
 `jprod_available`  | `Bool`        | true if the Jacobian-vector product `J * v` is available
 `jtprod_available` | `Bool`        | true if the transpose Jacobian-vector product `J' * v` is available
 `hprod_available`  | `Bool`        | true if the Hessian-vector product of the Lagrangian `H * v` is available
diff --git a/docs/src/reference.md b/docs/src/reference.md
index 45ffa7592..ef9a50ca1 100644
--- a/docs/src/reference.md
+++ b/docs/src/reference.md
@@ -1,17 +1,17 @@
 # Reference
-​
+
 ## Contents
-​
+
 ```@contents
 Pages = ["reference.md"]
 ```
-​
+
 ## Index
-​
+
 ```@index
 Pages = ["reference.md"]
 ```
-​
+
 ```@autodocs
 Modules = [NLPModels]
-```
\ No newline at end of file
+```
diff --git a/src/nlp/api.jl b/src/nlp/api.jl
index cd1efeb90..3b9933317 100644
--- a/src/nlp/api.jl
+++ b/src/nlp/api.jl
@@ -12,6 +12,7 @@ export jth_hprod, jth_hprod!, ghjvprod, ghjvprod!
 export hess_structure!, hess_structure, hess_coord!, hess_coord
 export hess, hprod, hprod!, hess_op, hess_op!
 export varscale, lagscale, conscale
+export jac_dense!, hess_dense!
 
 """
     f = obj(nlp, x)
@@ -56,26 +57,7 @@ end
 
 Evaluate ``c(x)``, the constraints at `x` in place.
 """
-function cons!(nlp::AbstractNLPModel, x::AbstractVector, cx::AbstractVector)
-  @lencheck nlp.meta.nvar x
-  @lencheck nlp.meta.ncon cx
-  increment!(nlp, :neval_cons)
-  if nlp.meta.nlin > 0
-    if nlp.meta.nnln == 0
-      cons_lin!(nlp, x, cx)
-    else
-      cons_lin!(nlp, x, view(cx, nlp.meta.lin))
-    end
-  end
-  if nlp.meta.nnln > 0
-    if nlp.meta.nlin == 0
-      cons_nln!(nlp, x, cx)
-    else
-      cons_nln!(nlp, x, view(cx, nlp.meta.nln))
-    end
-  end
-  return cx
-end
+function cons! end
 
 """
     c = cons_lin(nlp, x)
@@ -178,7 +160,7 @@ end
     (rows,cols) = jac_structure(nlp)
 
 Return the structure of the constraints Jacobian in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_structure(nlp::AbstractNLPModel)
   rows = Vector{Int}(undef, nlp.meta.nnzj)
@@ -190,44 +172,15 @@ end
     jac_structure!(nlp, rows, cols)
 
 Return the structure of the constraints Jacobian in sparse coordinate format in place.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
-function jac_structure!(
-  nlp::AbstractNLPModel,
-  rows::AbstractVector{T},
-  cols::AbstractVector{T},
-) where {T}
-  @lencheck nlp.meta.nnzj rows cols
-  lin_ind = 1:(nlp.meta.lin_nnzj)
-  if nlp.meta.nlin > 0
-    if nlp.meta.nnln == 0
-      jac_lin_structure!(nlp, rows, cols)
-    else
-      jac_lin_structure!(nlp, view(rows, lin_ind), view(cols, lin_ind))
-      for i in lin_ind
-        rows[i] += count(x < nlp.meta.lin[rows[i]] for x in nlp.meta.nln)
-      end
-    end
-  end
-  if nlp.meta.nnln > 0
-    if nlp.meta.nlin == 0
-      jac_nln_structure!(nlp, rows, cols)
-    else
-      nln_ind = (nlp.meta.lin_nnzj + 1):(nlp.meta.lin_nnzj + nlp.meta.nln_nnzj)
-      jac_nln_structure!(nlp, view(rows, nln_ind), view(cols, nln_ind))
-      for i in nln_ind
-        rows[i] += count(x < nlp.meta.nln[rows[i]] for x in nlp.meta.lin)
-      end
-    end
-  end
-  return rows, cols
-end
+function jac_structure! end
 
 """
     (rows,cols) = jac_lin_structure(nlp)
 
 Return the structure of the linear constraints Jacobian in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_lin_structure(nlp::AbstractNLPModel)
   rows = Vector{Int}(undef, nlp.meta.lin_nnzj)
@@ -239,7 +192,7 @@ end
     jac_lin_structure!(nlp, rows, cols)
 
 Return the structure of the linear constraints Jacobian in sparse coordinate format in place.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_lin_structure! end
 
@@ -247,7 +200,7 @@ function jac_lin_structure! end
     (rows,cols) = jac_nln_structure(nlp)
 
 Return the structure of the nonlinear constraints Jacobian in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_nln_structure(nlp::AbstractNLPModel)
   rows = Vector{Int}(undef, nlp.meta.nln_nnzj)
@@ -259,44 +212,22 @@ end
     jac_nln_structure!(nlp, rows, cols)
 
 Return the structure of the nonlinear constraints Jacobian in sparse coordinate format in place.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_nln_structure! end
 
 """
     vals = jac_coord!(nlp, x, vals)
 
-Evaluate ``J(x)``, the constraints Jacobian at `x` in sparse coordinate format, rewriting `vals`.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+Evaluate ``J(x)``, the constraints Jacobian at `x` in sparse coordinate format, overwriting `vals`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
-function jac_coord!(nlp::AbstractNLPModel, x::AbstractVector, vals::AbstractVector)
-  @lencheck nlp.meta.nvar x
-  @lencheck nlp.meta.nnzj vals
-  increment!(nlp, :neval_jac)
-  if nlp.meta.nlin > 0
-    if nlp.meta.nnln == 0
-      jac_lin_coord!(nlp, x, vals)
-    else
-      lin_ind = 1:(nlp.meta.lin_nnzj)
-      jac_lin_coord!(nlp, x, view(vals, lin_ind))
-    end
-  end
-  if nlp.meta.nnln > 0
-    if nlp.meta.nlin == 0
-      jac_nln_coord!(nlp, x, vals)
-    else
-      nln_ind = (nlp.meta.lin_nnzj + 1):(nlp.meta.lin_nnzj + nlp.meta.nln_nnzj)
-      jac_nln_coord!(nlp, x, view(vals, nln_ind))
-    end
-  end
-  return vals
-end
-
+function jac_coord! end
 """
     vals = jac_coord(nlp, x)
 
 Evaluate ``J(x)``, the constraints Jacobian at `x` in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_coord(nlp::AbstractNLPModel{T, S}, x::AbstractVector) where {T, S}
   @lencheck nlp.meta.nvar x
@@ -304,11 +235,19 @@ function jac_coord(nlp::AbstractNLPModel{T, S}, x::AbstractVector) where {T, S}
   return jac_coord!(nlp, x, vals)
 end
 
+"""
+    Jx = jac_dense!(nlp, x, Jx)
+
+Evaluate ``J(x)``, the constraints Jacobian at `x` in dense format, overwriting `Jx`.
+This function is only available when `nlp.meta.jac_available` is set to `true` and `nlp.meta.sparse_jacobian` is set to `false`.
+"""
+function jac_dense! end
+
 """
     Jx = jac(nlp, x)
 
 Evaluate ``J(x)``, the constraints Jacobian at `x` as a sparse matrix.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac(nlp::AbstractNLPModel, x::AbstractVector)
   @lencheck nlp.meta.nvar x
@@ -321,7 +260,7 @@ end
     vals = jac_lin_coord!(nlp, x, vals)
 
 Evaluate ``J(x)``, the linear constraints Jacobian at `x` in sparse coordinate format, overwriting `vals`.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_lin_coord! end
 
@@ -329,7 +268,7 @@ function jac_lin_coord! end
     vals = jac_lin_coord(nlp, x)
 
 Evaluate ``J(x)``, the linear constraints Jacobian at `x` in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_lin_coord(nlp::AbstractNLPModel{T, S}, x::AbstractVector) where {T, S}
   @lencheck nlp.meta.nvar x
@@ -341,7 +280,7 @@ end
     Jx = jac_lin(nlp, x)
 
 Evaluate ``J(x)``, the linear constraints Jacobian at `x` as a sparse matrix.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_lin(nlp::AbstractNLPModel, x::AbstractVector)
   @lencheck nlp.meta.nvar x
@@ -354,7 +293,7 @@ end
     vals = jac_nln_coord!(nlp, x, vals)
 
 Evaluate ``J(x)``, the nonlinear constraints Jacobian at `x` in sparse coordinate format, overwriting `vals`.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_nln_coord! end
 
@@ -362,7 +301,7 @@ function jac_nln_coord! end
     vals = jac_nln_coord(nlp, x)
 
 Evaluate ``J(x)``, the nonlinear constraints Jacobian at `x` in sparse coordinate format.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_nln_coord(nlp::AbstractNLPModel{T, S}, x::AbstractVector) where {T, S}
   @lencheck nlp.meta.nvar x
@@ -374,7 +313,7 @@ end
     Jx = jac_nln(nlp, x)
 
 Evaluate ``J(x)``, the nonlinear constraints Jacobian at `x` as a sparse matrix.
-This function is only available if `nlp.meta.jac_available` is set to `true`.
+This function is only available when both `nlp.meta.jac_available` and `nlp.meta.sparse_jacobian` are set to `true`.
 """
 function jac_nln(nlp::AbstractNLPModel, x::AbstractVector)
   @lencheck nlp.meta.nvar x
@@ -401,26 +340,7 @@ end
 Evaluate ``J(x)v``, the Jacobian-vector product at `x` in place.
 This function is only available if `nlp.meta.jprod_available` is set to `true`.
 """
-function jprod!(nlp::AbstractNLPModel, x::AbstractVector, v::AbstractVector, Jv::AbstractVector)
-  @lencheck nlp.meta.nvar x v
-  @lencheck nlp.meta.ncon Jv
-  increment!(nlp, :neval_jprod)
-  if nlp.meta.nlin > 0
-    if nlp.meta.nnln == 0
-      jprod_lin!(nlp, x, v, Jv)
-    else
-      jprod_lin!(nlp, x, v, view(Jv, nlp.meta.lin))
-    end
-  end
-  if nlp.meta.nnln > 0
-    if nlp.meta.nlin == 0
-      jprod_nln!(nlp, x, v, Jv)
-    else
-      jprod_nln!(nlp, x, v, view(Jv, nlp.meta.nln))
-    end
-  end
-  return Jv
-end
+function jprod! end
 
 """
     Jv = jprod!(nlp, rows, cols, vals, v, Jv)
@@ -545,27 +465,7 @@ Evaluate ``J(x)^Tv``, the transposed-Jacobian-vector product at `x` in place.
 If the problem has linear and nonlinear constraints, this function allocates.
 This function is only available if `nlp.meta.jtprod_available` is set to `true`.
 """
-function jtprod!(nlp::AbstractNLPModel, x::AbstractVector, v::AbstractVector, Jtv::AbstractVector)
-  @lencheck nlp.meta.nvar x Jtv
-  @lencheck nlp.meta.ncon v
-  increment!(nlp, :neval_jtprod)
-  if nlp.meta.nnln == 0
-    (nlp.meta.nlin > 0) && jtprod_lin!(nlp, x, v, Jtv)
-  elseif nlp.meta.nlin == 0
-    (nlp.meta.nnln > 0) && jtprod_nln!(nlp, x, v, Jtv)
-  elseif nlp.meta.nlin >= nlp.meta.nnln
-    jtprod_lin!(nlp, x, view(v, nlp.meta.lin), Jtv)
-    if nlp.meta.nnln > 0
-      Jtv .+= jtprod_nln(nlp, x, view(v, nlp.meta.nln))
-    end
-  else
-    jtprod_nln!(nlp, x, view(v, nlp.meta.nln), Jtv)
-    if nlp.meta.nlin > 0
-      Jtv .+= jtprod_lin(nlp, x, view(v, nlp.meta.lin))
-    end
-  end
-  return Jtv
-end
+function jtprod! end
 
 """
     Jtv = jtprod!(nlp, rows, cols, vals, v, Jtv)
@@ -953,6 +853,7 @@ end
 
 Evaluate the Hessian of j-th constraint at `x` in sparse coordinate format.
 Only the lower triangle is returned.
+This function is only available when `nlp.meta.sparse_hessian` is set to `true`.
 """
 function jth_hess_coord(nlp::AbstractNLPModel{T, S}, x::AbstractVector, j::Integer) where {T, S}
   @lencheck nlp.meta.nvar x
@@ -964,8 +865,9 @@ end
 """
     vals = jth_hess_coord!(nlp, x, j, vals)
 
-Evaluate the Hessian of j-th constraint at `x` in sparse coordinate format, with `vals` of
-length `nlp.meta.nnzh`, in place. Only the lower triangle is returned.
+Evaluate the Hessian of j-th constraint at `x` in sparse coordinate format, with `vals` of length `nlp.meta.nnzh`, in place.
+Only the lower triangle is returned.
+This function is only available when `nlp.meta.sparse_hessian` is set to `true`.
 """
 function jth_hess_coord! end
 
@@ -975,6 +877,7 @@ function jth_hess_coord! end
 Evaluate the Hessian of j-th constraint at `x` as a sparse matrix with
 the same sparsity pattern as the Lagrangian Hessian.
 A `Symmetric` object wrapping the lower triangle is returned.
+This function is only available when `nlp.meta.sparse_hessian` is set to `true`.
 """
 function jth_hess(nlp::AbstractNLPModel, x::AbstractVector, j::Integer)
   @lencheck nlp.meta.nvar x
@@ -1036,7 +939,7 @@ function ghjvprod! end
     (rows,cols) = hess_structure(nlp)
 
 Return the structure of the Lagrangian Hessian in sparse coordinate format.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_structure(nlp::AbstractNLPModel)
   rows = Vector{Int}(undef, nlp.meta.nnzh)
@@ -1048,7 +951,7 @@ end
     hess_structure!(nlp, rows, cols)
 
 Return the structure of the Lagrangian Hessian in sparse coordinate format in place.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_structure! end
 
@@ -1059,7 +962,7 @@ Evaluate the objective Hessian at `x` in sparse coordinate format,
 with objective function scaled by `obj_weight`, i.e.,
 $(OBJECTIVE_HESSIAN), overwriting `vals`.
 Only the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_coord!(
   nlp::AbstractNLPModel{T, S},
@@ -1080,7 +983,7 @@ Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
 with objective function scaled by `obj_weight`, i.e.,
 $(LAGRANGIAN_HESSIAN), overwriting `vals`.
 Only the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_coord! end
 
@@ -1091,7 +994,7 @@ Evaluate the objective Hessian at `x` in sparse coordinate format,
 with objective function scaled by `obj_weight`, i.e.,
 $(OBJECTIVE_HESSIAN).
 Only the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_coord(
   nlp::AbstractNLPModel{T, S},
@@ -1110,7 +1013,7 @@ Evaluate the Lagrangian Hessian at `(x,y)` in sparse coordinate format,
 with objective function scaled by `obj_weight`, i.e.,
 $(LAGRANGIAN_HESSIAN).
 Only the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess_coord(
   nlp::AbstractNLPModel{T, S},
@@ -1131,7 +1034,7 @@ Evaluate the objective Hessian at `x` as a sparse matrix,
 with objective function scaled by `obj_weight`, i.e.,
 $(OBJECTIVE_HESSIAN).
 A `Symmetric` object wrapping the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess(
   nlp::AbstractNLPModel{T, S},
@@ -1144,6 +1047,17 @@ function hess(
   Symmetric(sparse(rows, cols, vals, nlp.meta.nvar, nlp.meta.nvar), :L)
 end
 
+"""
+    Hx = hess_dense!(nlp, x, Hx; obj_weight=1.0)
+    Hx = hess_dense!(nlp, x, y, Hx; obj_weight=1.0)
+
+The first method evaluates ``H(x)``, the Hessian of the objective at `x` in dense format, overwriting `Hx`.
+The second method evaluates ``H(x,y)``, the Hessian of the Lagrangian at `(x,y)` in dense format, overwriting `Hx`.
+Only the lower triangular part of `Hx` needs to be filled.
+This function is only available when `nlp.meta.hess_available` is set to `true` and `nlp.meta.sparse_hessian` is set to `false`.
+"""
+function hess_dense! end
+
 """
     Hx = hess(nlp, x, y; obj_weight=1.0)
 
@@ -1151,7 +1065,7 @@ Evaluate the Lagrangian Hessian at `(x,y)` as a sparse matrix,
 with objective function scaled by `obj_weight`, i.e.,
 $(LAGRANGIAN_HESSIAN).
 A `Symmetric` object wrapping the lower triangle is returned.
-This function is only available if `nlp.meta.hess_available` is set to `true`.
+This function is only available when both `nlp.meta.hess_available` and `nlp.meta.sparse_hessian` are set to `true`.
 """
 function hess(
   nlp::AbstractNLPModel{T, S},
diff --git a/src/nlp/meta.jl b/src/nlp/meta.jl
index 884d0922a..491c49dfa 100644
--- a/src/nlp/meta.jl
+++ b/src/nlp/meta.jl
@@ -43,17 +43,19 @@ The following keyword arguments are accepted:
 - `lcon`: vector of constraint lower bounds
 - `ucon`: vector of constraint upper bounds
 - `nnzo`: number of nonzeros in the gradient
-- `nnzj`: number of elements needed to store the nonzeros in the sparse Jacobian
-- `lin_nnzj`: number of elements needed to store the nonzeros in the sparse Jacobian of linear constraints
-- `nln_nnzj`: number of elements needed to store the nonzeros in the sparse Jacobian of nonlinear constraints
-- `nnzh`: number of elements needed to store the nonzeros in the sparse Hessian
+- `nnzj`: number of elements needed to store the nonzeros in the Jacobian
+- `lin_nnzj`: number of elements needed to store the nonzeros in the Jacobian of linear constraints
+- `nln_nnzj`: number of elements needed to store the nonzeros in the Jacobian of nonlinear constraints
+- `nnzh`: number of elements needed to store the nonzeros in the Hessian of the Lagrangian
 - `lin`: indices of linear constraints
 - `minimize`: true if optimize == minimize
 - `islp`: true if the problem is a linear program
 - `name`: problem name
+- `sparse_jacobian`: indicates whether the Jacobian of the constraints is sparse
+- `sparse_hessian`: indicates whether the Hessian of the Lagrangian is sparse
 - `grad_available`: indicates whether the gradient of the objective is available
-- `jac_available`: indicates whether the sparse Jacobian of the constraints is available
-- `hess_available`: indicates whether the sparse Hessian of the Lagrangian is available
+- `jac_available`: indicates whether the Jacobian of the constraints is available
+- `hess_available`: indicates whether the Hessian of the Lagrangian is available
 - `jprod_available`: indicates whether the Jacobian-vector product `J * v` is available
 - `jtprod_available`: indicates whether the transpose Jacobian-vector product `J' * v` is available
 - `hprod_available`: indicates whether the Hessian-vector product of the Lagrangian `H * v` is available
@@ -121,6 +123,9 @@ struct NLPModelMeta{T, S} <: AbstractNLPModelMeta{T, S}
   islp::Bool
   name::String
 
+  sparse_jacobian::Bool
+  sparse_hessian::Bool
+
   grad_available::Bool
   jac_available::Bool
   hess_available::Bool
@@ -145,11 +150,13 @@ function NLPModelMeta{T, S}(
   nnzj = nvar * ncon,
   lin_nnzj = 0,
   nln_nnzj = nnzj - lin_nnzj,
-  nnzh = nvar * (nvar + 1) / 2,
+  nnzh = nvar * (nvar + 1) ÷ 2,
   lin = Int[],
   minimize::Bool = true,
   islp::Bool = false,
   name = "Generic",
+  sparse_jacobian::Bool = true,
+  sparse_hessian::Bool = true,
   grad_available::Bool = true,
   jac_available::Bool = (ncon > 0),
   hess_available::Bool = true,
@@ -157,7 +164,7 @@ function NLPModelMeta{T, S}(
   jtprod_available::Bool = (ncon > 0),
   hprod_available::Bool = true,
 ) where {T, S}
-  if (nvar < 1) || (ncon < 0)
+  if (nvar < 1) || (ncon < 0) || (nnzj < 0) || (nnzh < 0)
     error("Nonsensical dimensions")
   end
 
@@ -189,9 +196,6 @@ function NLPModelMeta{T, S}(
     jinf = Int[]
   end
 
-  nnzj = max(0, nnzj)
-  nnzh = max(0, nnzh)
-
   nln = setdiff(1:ncon, lin)
   nlin = length(lin)
   nnln = length(nln)
@@ -232,6 +236,8 @@ function NLPModelMeta{T, S}(
     minimize,
     islp,
     name,
+    sparse_jacobian,
+    sparse_hessian,
     grad_available,
     jac_available,
     hess_available,
@@ -266,6 +272,8 @@ function NLPModelMeta(
   minimize::Bool = meta.minimize,
   islp::Bool = meta.islp,
   name = meta.name,
+  sparse_jacobian::Bool = true,
+  sparse_hessian::Bool = true,
   grad_available::Bool = meta.grad_available,
   jac_available::Bool = meta.jac_available,
   hess_available::Bool = meta.hess_available,
@@ -294,6 +302,8 @@ function NLPModelMeta(
     minimize = minimize,
     islp = islp,
     name = name,
+    sparse_jacobian = sparse_jacobian,
+    sparse_hessian = sparse_hessian,
     grad_available = grad_available,
     jac_available = jac_available,
     hess_available = hess_available,
diff --git a/test/nlp/simple-model.jl b/test/nlp/simple-model.jl
index 7f61c6260..745bd6c75 100644
--- a/test/nlp/simple-model.jl
+++ b/test/nlp/simple-model.jl
@@ -94,11 +94,20 @@ function NLPModels.hprod!(
   return Hv
 end
 
+function NLPModels.cons!(nlp::SimpleNLPModel, x::AbstractVector, cx::AbstractVector)
+  @lencheck 2 x
+  @lencheck 2 cx
+  increment!(nlp, :neval_cons)
+  cx[1] = x[1] - 2 * x[2] + 1
+  cx[2] = -x[1]^2 / 4 - x[2]^2 + 1
+  return cx
+end
+
 function NLPModels.cons_nln!(nlp::SimpleNLPModel, x::AbstractVector, cx::AbstractVector)
   @lencheck 2 x
   @lencheck 1 cx
   increment!(nlp, :neval_cons_nln)
-  cx .= [-x[1]^2 / 4 - x[2]^2 + 1]
+  cx[1] = -x[1]^2 / 4 - x[2]^2 + 1
   return cx
 end
 
@@ -106,10 +115,21 @@ function NLPModels.cons_lin!(nlp::SimpleNLPModel, x::AbstractVector, cx::Abstrac
   @lencheck 2 x
   @lencheck 1 cx
   increment!(nlp, :neval_cons_lin)
-  cx .= [x[1] - 2 * x[2] + 1]
+  cx[1] = x[1] - 2 * x[2] + 1
   return cx
 end
 
+function NLPModels.jac_structure!(
+  nlp::SimpleNLPModel,
+  rows::AbstractVector{Int},
+  cols::AbstractVector{Int},
+)
+  @lencheck 4 rows cols
+  rows .= [1, 1, 2, 2]
+  cols .= [1, 2, 1, 2]
+  return rows, cols
+end
+
 function NLPModels.jac_nln_structure!(
   nlp::SimpleNLPModel,
   rows::AbstractVector{Int},
@@ -132,20 +152,48 @@ function NLPModels.jac_lin_structure!(
   return rows, cols
 end
 
+function NLPModels.jac_coord!(nlp::SimpleNLPModel, x::AbstractVector, vals::AbstractVector)
+  @lencheck 2 x
+  @lencheck 4 vals
+  increment!(nlp, :neval_jac)
+  vals[1] = 1
+  vals[2] = -2
+  vals[3] = -x[1] / 2
+  vals[4] = -2 * x[2]
+  return vals
+end
+
 function NLPModels.jac_nln_coord!(nlp::SimpleNLPModel, x::AbstractVector, vals::AbstractVector)
-  @lencheck 2 x vals
+  @lencheck 2 x
+  @lencheck 2 vals
   increment!(nlp, :neval_jac_nln)
-  vals .= [-x[1] / 2, -2 * x[2]]
+  vals[1] = -x[1] / 2
+  vals[2] = -2 * x[2]
   return vals
 end
 
 function NLPModels.jac_lin_coord!(nlp::SimpleNLPModel, x::AbstractVector, vals::AbstractVector)
   @lencheck 2 x vals
   increment!(nlp, :neval_jac_lin)
-  vals .= [1, -2]
+  vals[1] = 1
+  vals[2] = -2
   return vals
 end
 
+function NLPModels.jprod!(
+  nlp::SimpleNLPModel,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jv::AbstractVector,
+)
+  @lencheck 2 x v
+  @lencheck 2 Jv
+  increment!(nlp, :neval_jprod)
+  Jv[1] = v[1] - 2 * v[2]
+  Jv[2] = -x[1] * v[1] / 2 - 2 * x[2] * v[2]
+  return Jv
+end
+
 function NLPModels.jprod_nln!(
   nlp::SimpleNLPModel,
   x::AbstractVector,
@@ -155,7 +203,7 @@ function NLPModels.jprod_nln!(
   @lencheck 2 x v
   @lencheck 1 Jv
   increment!(nlp, :neval_jprod_nln)
-  Jv .= [-x[1] * v[1] / 2 - 2 * x[2] * v[2]]
+  Jv[1] = -x[1] * v[1] / 2 - 2 * x[2] * v[2]
   return Jv
 end
 
@@ -168,10 +216,24 @@ function NLPModels.jprod_lin!(
   @lencheck 2 x v
   @lencheck 1 Jv
   increment!(nlp, :neval_jprod_lin)
-  Jv .= [v[1] - 2 * v[2]]
+  Jv[1] = v[1] - 2 * v[2]
   return Jv
 end
 
+function NLPModels.jtprod!(
+  nlp::SimpleNLPModel,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jtv::AbstractVector,
+)
+  @lencheck 2 x Jtv
+  @lencheck 2 v
+  increment!(nlp, :neval_jtprod)
+  Jtv[1] = v[1] - x[1] * v[2] / 2
+  Jtv[2] = -2 * v[1] -2 * x[2] * v[2]
+  return Jtv
+end
+
 function NLPModels.jtprod_nln!(
   nlp::SimpleNLPModel,
   x::AbstractVector,
@@ -181,7 +243,8 @@ function NLPModels.jtprod_nln!(
   @lencheck 2 x Jtv
   @lencheck 1 v
   increment!(nlp, :neval_jtprod_nln)
-  Jtv .= [-x[1] * v[1] / 2; -2 * x[2] * v[1]]
+  Jtv[1] = -x[1] * v[1] / 2
+  Jtv[2] = -2 * x[2] * v[1]
   return Jtv
 end
 
@@ -194,7 +257,8 @@ function NLPModels.jtprod_lin!(
   @lencheck 2 x Jtv
   @lencheck 1 v
   increment!(nlp, :neval_jtprod_lin)
-  Jtv .= [v[1]; -2 * v[1]]
+  Jtv[1] = v[1]
+  Jtv[2] = -2 * v[1]
   return Jtv
 end
 
diff --git a/test/nls/simple-model.jl b/test/nls/simple-model.jl
index 8f1fe50c8..56d9672bd 100644
--- a/test/nls/simple-model.jl
+++ b/test/nls/simple-model.jl
@@ -135,6 +135,14 @@ function NLPModels.hprod_residual!(
   return Hiv
 end
 
+function NLPModels.cons!(nls::SimpleNLSModel, x::AbstractVector, cx::AbstractVector)
+  @lencheck 2 x
+  @lencheck 3 cx
+  increment!(nls, :neval_cons)
+  cx .= [x[1] + x[2]^2; x[1]^2 + x[2]; x[1]^2 + x[2]^2 - 1]
+  return cx
+end
+
 function NLPModels.cons_nln!(nls::SimpleNLSModel, x::AbstractVector, cx::AbstractVector)
   @lencheck 2 x
   @lencheck 3 cx
@@ -143,6 +151,17 @@ function NLPModels.cons_nln!(nls::SimpleNLSModel, x::AbstractVector, cx::Abstrac
   return cx
 end
 
+function NLPModels.jac_structure!(
+  nls::SimpleNLSModel,
+  rows::AbstractVector{<:Integer},
+  cols::AbstractVector{<:Integer},
+)
+  @lencheck 6 rows cols
+  rows .= [1, 1, 2, 2, 3, 3]
+  cols .= [1, 2, 1, 2, 1, 2]
+  return rows, cols
+end
+
 function NLPModels.jac_nln_structure!(
   nls::SimpleNLSModel,
   rows::AbstractVector{<:Integer},
@@ -154,6 +173,14 @@ function NLPModels.jac_nln_structure!(
   return rows, cols
 end
 
+function NLPModels.jac_coord!(nls::SimpleNLSModel, x::AbstractVector, vals::AbstractVector)
+  @lencheck 2 x
+  @lencheck 6 vals
+  increment!(nls, :neval_jac)
+  vals .= [1, 2x[2], 2x[1], 1, 2x[1], 2x[2]]
+  return vals
+end
+
 function NLPModels.jac_nln_coord!(nls::SimpleNLSModel, x::AbstractVector, vals::AbstractVector)
   @lencheck 2 x
   @lencheck 6 vals
@@ -162,6 +189,19 @@ function NLPModels.jac_nln_coord!(nls::SimpleNLSModel, x::AbstractVector, vals::
   return vals
 end
 
+function NLPModels.jprod!(
+  nls::SimpleNLSModel,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jv::AbstractVector,
+)
+  @lencheck 2 x v
+  @lencheck 3 Jv
+  increment!(nls, :neval_jprod)
+  Jv .= [v[1] + 2x[2] * v[2]; 2x[1] * v[1] + v[2]; 2x[1] * v[1] + 2x[2] * v[2]]
+  return Jv
+end
+
 function NLPModels.jprod_nln!(
   nls::SimpleNLSModel,
   x::AbstractVector,
@@ -175,6 +215,19 @@ function NLPModels.jprod_nln!(
   return Jv
 end
 
+function NLPModels.jtprod!(
+  nls::SimpleNLSModel,
+  x::AbstractVector,
+  v::AbstractVector,
+  Jtv::AbstractVector,
+)
+  @lencheck 2 x Jtv
+  @lencheck 3 v
+  increment!(nls, :neval_jtprod)
+  Jtv .= [v[1] + 2x[1] * (v[2] + v[3]); v[2] + 2x[2] * (v[1] + v[3])]
+  return Jtv
+end
+
 function NLPModels.jtprod_nln!(
   nls::SimpleNLSModel,
   x::AbstractVector,