diff --git a/README.md b/README.md
index 3557944b..a6781c55 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 c7424fb8..cd91de80 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 b5e42634..8562e197 100644
--- a/docs/src/guidelines.md
+++ b/docs/src/guidelines.md
@@ -60,21 +60,24 @@ 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)`
+  - `jac_lin_structure!(nlp, jrows, jcols)` (sparse Jacobian)
+  - `jac_nln_structure!(nlp, jrows, jcols)` (sparse Jacobian)
+  - `jac_lin_coord!(nlp, x, jvals)` (sparse Jacobian)
+  - `jac_nln_coord!(nlp, x, jvals)` (sparse Jacobian)
+  - `jac_dense!(nlp, x, Jx)` (dense Jacobian)
   - `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)`
+  - `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)`
 
 The linear constraints are specified at the initialization of the `NLPModelMeta` using the keyword arguement `lin`.
@@ -82,6 +85,12 @@ The indices of linear and nonlinear constraints are respectively available in `n
 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.
+This is specified at the initialization of [`NLPModelMeta`](@ref) through the keyword arguments `sparse_jacobian` and `sparse_hessian`.
+In the dense case, linear and nonlinear constraints are handled together.
+Only an in-place API is available for dense Jacobians and Hessians (`jac_dense!` and `hess_dense!`).
+
 ## [Availability of the API](@id availability-api)
 
 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
diff --git a/docs/src/index.md b/docs/src/index.md
index c5c4f123..3a8758cc 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/src/nlp/api.jl b/src/nlp/api.jl
index cd1efeb9..ca15e375 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)
@@ -178,7 +179,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,7 +191,7 @@ 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,
@@ -227,7 +228,7 @@ 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 +240,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 +248,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,15 +260,15 @@ 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
@@ -296,7 +297,7 @@ 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 +305,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 +330,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 +338,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 +350,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 +363,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 +371,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 +383,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
@@ -953,6 +962,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 +974,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 +986,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 +1048,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 +1060,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 +1071,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 +1092,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 +1103,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 +1122,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 +1143,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 +1156,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 +1174,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 884d0922..491c49df 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,