diff --git a/docs/src/submodules/Utilities/overview.md b/docs/src/submodules/Utilities/overview.md
index 084e746c82..4192bfd9a6 100644
--- a/docs/src/submodules/Utilities/overview.md
+++ b/docs/src/submodules/Utilities/overview.md
@@ -11,6 +11,7 @@ DocTestFilters = [r"MathOptInterface|MOI"]
 
 The Utilities submodule provides a variety of functionality for managing
 `MOI.ModelLike` objects.
+
 ## Utilities.Model
 
 [`Utilities.Model`](@ref) provides an implementation of a [`ModelLike`](@ref)
@@ -382,6 +383,93 @@ See [`Utilities.supports_allocate_load`](@ref) for more details.
 !!! warning
     The Allocate-Load API should **not** be used outside [`copy_to`](@ref).
 
+## Utilities.MatrixOfConstraints
+
+The constraints of [`Utilities.Model`](@ref) are stored as a vector of tuples
+of function and set in a `Utilities.VectorOfConstraints`. Other representations
+can be used by parametrizing the type [`Utilities.GenericModel`](@ref)
+(resp. [`Utilities.GenericOptimizer`](@ref)). For instance, if all
+non-`SingleVariable` constraints are affine, the coefficients of all the
+constraints can be stored in a single sparse matrix using
+[`Utilities.MatrixOfConstraints`](@ref). The constraints storage can even be
+customized up to a point where it exactly matches the storage of the solver of
+interest, in which case [`copy_to`](@ref) can be implemented for the solver by
+calling [`copy_to`](@ref) to this custom model.
+
+For instance, [Clp](https://github.com/jump-dev/Clp.jl) defines the following
+model
+```julia
+MOI.Utilities.@product_of_scalar_sets(LP, MOI.EqualTo{T}, MOI.LessThan{T}, MOI.GreaterThan{T})
+const Model = MOI.Utilities.GenericModel{
+    Float64,
+    MOI.Utilities.MatrixOfConstraints{
+        Float64,
+        MOI.Utilities.MutableSparseMatrixCSC{Float64,Cint,MOI.Utilities.ZeroBasedIndexing},
+        MOI.Utilities.Box{Float64},
+        LP{Float64},
+    },
+}
+```
+
+The [`copy_to`](@ref) operation can now be implemented as follows (assuming that
+the `Model` definition above is in the `Clp` module so that it can be referred
+to as `Model`, to be distinguished with [`Utilities.Model`](@ref)):
+```julia
+function _copy_to(
+    dest::Optimizer,
+    src::Model
+)
+    @assert MOI.is_empty(dest)
+    A = src.constraints.coefficients
+    row_bounds = src.constraints.constants
+    Clp_loadProblem(
+        dest,
+        A.n,
+        A.m,
+        A.colptr,
+        A.rowval,
+        A.nzval,
+        src.lower_bound,
+        src.upper_bound,
+        # (...) objective vector (omitted),
+        row_bounds.lower,
+        row_bounds.upper,
+    )
+    # Set objective sense and constant (omitted)
+    return
+end
+
+function MOI.copy_to(
+    dest::Optimizer,
+    src::Model;
+    copy_names::Bool = false
+)
+    _copy_to(dest, src)
+    return MOI.Utilities.identity_index_map(src)
+end
+
+function MOI.copy_to(
+    dest::Optimizer,
+    src::MOI.Utilities.UniversalFallback{Model};
+    copy_names::Bool = false
+)
+    # Copy attributes from `src` to `dest` and error in case any unsupported
+    # constraints or attributes are set in `UniversalFallback`.
+    return MOI.copy_to(dest, src.model)
+end
+
+function MOI.copy_to(
+    dest::Optimizer,
+    src::MOI.ModelLike;
+    copy_names::Bool = false
+)
+    model = Model()
+    index_map = MOI.copy_to(model, src)
+    _copy_to(dest, model)
+    return index_map
+end
+```
+
 ## Fallbacks
 
 The value of some attributes can be inferred from the value of other
diff --git a/docs/src/submodules/Utilities/reference.md b/docs/src/submodules/Utilities/reference.md
index fbb55a27d3..a70a845cd5 100644
--- a/docs/src/submodules/Utilities/reference.md
+++ b/docs/src/submodules/Utilities/reference.md
@@ -79,6 +79,47 @@ Utilities.load
 Utilities.load_constraint
 ```
 
+## Utilities.MatrixOfConstraints
+
+```@docs
+Utilities.MatrixOfConstraints
+Utilities.rows
+```
+
+### Constants
+
+```@docs
+Utilities.Box
+Utilities.load_constants
+```
+
+### Mutable sparse matrix
+
+```@docs
+Utilities.MutableSparseMatrixCSC
+Utilities.AbstractIndexing
+Utilities.ZeroBasedIndexing
+Utilities.OneBasedIndexing
+Utilities.add_column
+Utilities.allocate_terms
+Utilities.set_number_of_rows
+Utilities.load_terms
+Utilities.final_touch
+```
+
+### Product of sets
+
+```@docs
+Utilities.ProductOfSets
+Utilities.set_index
+Utilities.set_types
+Utilities.add_set
+Utilities.indices
+Utilities.MixOfScalarSets
+Utilities.OrderedProductOfSets
+Utilities.OrderedProductOfScalarSets
+```
+
 ## Fallbacks
 
 ```@docs
diff --git a/src/Bridges/bridge_optimizer.jl b/src/Bridges/bridge_optimizer.jl
index d7cbfcb26f..08d1838bc3 100644
--- a/src/Bridges/bridge_optimizer.jl
+++ b/src/Bridges/bridge_optimizer.jl
@@ -395,6 +395,9 @@ function MOI.supports_incremental_interface(
 )
     return MOI.supports_incremental_interface(b.model, copy_names)
 end
+function MOIU.final_touch(uf::AbstractBridgeOptimizer, index_map)
+    return MOIU.final_touch(uf.model, index_map)
+end
 
 # References
 function MOI.is_valid(b::AbstractBridgeOptimizer, vi::MOI.VariableIndex)
diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl
index b689f2f0d1..69e47a3e38 100644
--- a/src/Utilities/Utilities.jl
+++ b/src/Utilities/Utilities.jl
@@ -60,6 +60,9 @@ include("variables.jl")
 include("vector_of_constraints.jl")
 include("struct_of_constraints.jl")
 include("model.jl")
+include("sparse_matrix.jl")
+include("product_of_sets.jl")
+include("matrix_of_constraints.jl")
 include("parser.jl")
 include("mockoptimizer.jl")
 include("cachingoptimizer.jl")
diff --git a/src/Utilities/cachingoptimizer.jl b/src/Utilities/cachingoptimizer.jl
index c292586311..0966db19c6 100644
--- a/src/Utilities/cachingoptimizer.jl
+++ b/src/Utilities/cachingoptimizer.jl
@@ -232,6 +232,27 @@ end
 
 _standardize(d::IndexMap) = d
 
+function pass_nonvariable_constraints(
+    dest::CachingOptimizer,
+    src::MOI.ModelLike,
+    idxmap::IndexMap,
+    constraint_types,
+    pass_cons;
+    filter_constraints::Union{Nothing,Function} = nothing,
+)
+    dest.state == ATTACHED_OPTIMIZER && reset_optimizer(dest)
+    return pass_nonvariable_constraints(
+        dest.model_cache,
+        src,
+        idxmap,
+        constraint_types,
+        pass_cons;
+        filter_constraints = filter_constraints,
+    )
+end
+function final_touch(m::CachingOptimizer, index_map)
+    return final_touch(m.model_cache, index_map)
+end
 function MOI.copy_to(m::CachingOptimizer, src::MOI.ModelLike; kws...)
     if m.state == ATTACHED_OPTIMIZER
         reset_optimizer(m)
diff --git a/src/Utilities/copy.jl b/src/Utilities/copy.jl
index 9d052820f9..64f4b654fb 100644
--- a/src/Utilities/copy.jl
+++ b/src/Utilities/copy.jl
@@ -569,6 +569,15 @@ function try_constrain_variables_on_creation(
     single_variable_not_added
 end
 
+"""
+    function final_touch(model::MOI.ModelLike, idxmap) end
+
+This is called at the end of [`default_copy_to`](@ref) to inform the model that
+the copy is finished. This allows `model` to perform thats that should be done
+only once all the model information is gathered.
+"""
+function final_touch(::MOI.ModelLike, idxmap) end
+
 """
     default_copy_to(
         dest::MOI.ModelLike,
@@ -649,6 +658,8 @@ function default_copy_to(
         filter_constraints = filter_constraints,
     )
 
+    final_touch(dest, idxmap)
+
     return idxmap
 end
 
diff --git a/src/Utilities/matrix_of_constraints.jl b/src/Utilities/matrix_of_constraints.jl
new file mode 100644
index 0000000000..622a0eafa7
--- /dev/null
+++ b/src/Utilities/matrix_of_constraints.jl
@@ -0,0 +1,316 @@
+# Constants: stored in a `Vector` or `Box` or any other type implementing:
+# `empty!`, `resize!` and `load_constants`.
+
+"""
+    load_constants(constants, offset, func_or_set)
+
+This function loads the constants of `func_or_set` in `constants` at an offset
+of `offset`. Where `offset` is the sum of the dimensions of the constraints
+already loaded. The storage should be preallocated with `resize!` before calling
+this function.
+
+This function should be implemented to be usable as storage of constants for
+[`MatrixOfConstraints`](@ref).
+
+# The constants are loaded in three steps:
+# 1) `Base.empty!` is called.
+# 2) `Base.resize!` is called with the sum of the dimensions of all constraints.
+# 3) `MOI.Utilities.load_constants` is called for each function for vector
+#    constraint or set for scalar constraint.
+"""
+function load_constants end
+
+function load_constants(
+    b::Vector{T},
+    offset,
+    func::MOI.VectorAffineFunction{T},
+) where {T}
+    copyto!(b, offset + 1, func.constants)
+    return
+end
+
+"""
+    struct Box{T}
+        lower::Vector{T}
+        upper::Vector{T}
+    end
+
+Stores the constants of scalar constraints with the lower bound of the set in
+`lower` and the upper bound in `upper`.
+"""
+struct Box{T}
+    lower::Vector{T}
+    upper::Vector{T}
+end
+Box{T}() where {T} = Box{T}(T[], T[])
+Base.:(==)(a::Box, b::Box) = a.lower == b.lower && a.upper == b.upper
+function Base.empty!(b::Box)
+    empty!(b.lower)
+    empty!(b.upper)
+    return b
+end
+
+function Base.resize!(b::Box, n)
+    resize!(b.lower, n)
+    resize!(b.upper, n)
+    return
+end
+
+function load_constants(
+    b::Box{T},
+    offset,
+    set::SUPPORTED_VARIABLE_SCALAR_SETS{T},
+) where {T}
+    flag = single_variable_flag(typeof(set))
+    b.lower[offset+1] = if iszero(flag & LOWER_BOUND_MASK)
+        typemin(T)
+    else
+        extract_lower_bound(set)
+    end
+    b.upper[offset+1] = if iszero(flag & UPPER_BOUND_MASK)
+        typemax(T)
+    else
+        extract_upper_bound(set)
+    end
+    return
+end
+
+"""
+    mutable struct MatrixOfConstraints{T,AT,BT,ST} <: MOI.ModelLike
+        coefficients::AT
+        constants::BT
+        sets::ST
+        are_indices_mapped::BitSet
+        caches::Union{Nothing, Vector}
+        function MatrixOfConstraints{T,AT,BT,ST}() where {T,AT,BT,ST}
+            return new{T,AT,BT,ST}(AT(), BT(), ST(), BitSet(), nothing)
+        end
+    end
+
+Represents affine constraints in a matrix form where the linear coefficients
+of the functions are stored in the `coefficients` field and the constants of the
+functions or sets are stored in the `sets` field. Additional information
+about the sets are stored in the `sets` field.
+
+This model can only be used as the `constraints` field of a
+`MOI.Utilities.AbstractModel`. When the constraints are added,
+they are stored in the `caches` field. They are only loaded in
+the `coefficients` and `constants` fields once `MOI.Utilities.final_touch`
+is called. For this reason, this should not be used with incremental
+building of the model but with a `MOI.copy_to` instead.
+
+The constraints can be added in two different ways:
+1) With `add_constraint` in which case a canonicalized copy
+   of the function is stored in `caches`.
+2) With `pass_nonvariable_constraints` in which case the functions and sets are
+   stored themselves in `caches` without mapping the variable indices.
+   The corresponding index in `caches` is added in `are_indices_mapped`.
+   This allows to avoid doing a copy of the function in case
+   the getter of `CanonicalConstraintFunction` does not make a copy
+   for the source model, e.g., this is the case of `VectorOfConstraints`.
+
+We illustrate this with an example. Suppose a model is copied from
+a `src::MOI.Utilities.Model` to a bridged model with a `MatrixOfConstraints`.
+For all the types that are not bridged, the constraints will be copied
+with `pass_nonvariable_constraints` hence the functions stored in
+`caches` are exactly the same as the ones stored in `src`.
+This is ok since this is only during the `copy_to` operation during which `src`
+cannot be modified.
+On the other hand, for the types that are bridged, the functions added
+may contain duplicates even if the functions did not contain duplicates in
+`src` so duplicates are removed with `MOI.Utilities.canonical`.
+"""
+mutable struct MatrixOfConstraints{T,AT,BT,ST} <: MOI.ModelLike
+    coefficients::AT
+    constants::BT
+    sets::ST
+    are_indices_mapped::BitSet
+    caches::Union{Nothing,Vector}
+    function MatrixOfConstraints{T,AT,BT,ST}() where {T,AT,BT,ST}
+        return new{T,AT,BT,ST}(AT(), BT(), ST(), BitSet(), nothing)
+    end
+end
+
+MOI.is_empty(v::MatrixOfConstraints) = MOI.is_empty(v.sets)
+function MOI.empty!(v::MatrixOfConstraints{T}) where {T}
+    MOI.empty!(v.coefficients)
+    empty!(v.constants)
+    MOI.empty!(v.sets)
+    empty!(v.are_indices_mapped)
+    v.caches =
+        [Tuple{_affine_function_type(T, S),S}[] for S in set_types(v.sets)]
+    return
+end
+
+"""
+    rows(model::MatrixOfConstraints, ci::MOI.ConstraintIndex)
+
+Return the rows corresponding to the constraint of index `ci`. If it is a
+vector constraint, this is a `UnitRange`, otherwise, this is an integer.
+"""
+function rows(model::MatrixOfConstraints, ci::MOI.ConstraintIndex)
+    return indices(model.sets, ci)
+end
+
+function _affine_function_type(
+    ::Type{T},
+    ::Type{<:MOI.AbstractScalarSet},
+) where {T}
+    return MOI.ScalarAffineFunction{T}
+end
+
+function _affine_function_type(
+    ::Type{T},
+    ::Type{<:MOI.AbstractVectorSet},
+) where {T}
+    return MOI.VectorAffineFunction{T}
+end
+
+function MOI.supports_constraint(
+    v::MatrixOfConstraints{T},
+    ::Type{F},
+    ::Type{S},
+) where {T,F<:MOI.AbstractFunction,S<:MOI.AbstractSet}
+    return F == _affine_function_type(T, S) && set_index(v.sets, S) !== nothing
+end
+
+function MOI.is_valid(v::MatrixOfConstraints, ci::MOI.ConstraintIndex)
+    return MOI.is_valid(v.sets, ci)
+end
+
+function MOI.get(
+    v::MatrixOfConstraints,
+    attr::Union{
+        MOI.ListOfConstraintTypesPresent,
+        MOI.NumberOfConstraints,
+        MOI.ListOfConstraintIndices,
+    },
+)
+    return MOI.get(v.sets, attr)
+end
+
+_add_set(sets, i, func::MOI.AbstractScalarFunction) = add_set(sets, i)
+function _add_set(sets, i, func::MOI.AbstractVectorFunction)
+    return add_set(sets, i, MOI.output_dimension(func))
+end
+function _add_constraint(model::MatrixOfConstraints, i, index_map, func, set)
+    allocate_terms(model.coefficients, index_map, func)
+    # Without this type annotation, the compiler is unable to know the type
+    # of `caches[i]` so this is slower and produce an allocation.
+    push!(model.caches[i]::Vector{Tuple{typeof(func),typeof(set)}}, (func, set))
+    return MOI.ConstraintIndex{typeof(func),typeof(set)}(
+        _add_set(model.sets, i, func),
+    )
+end
+
+struct IdentityMap <: AbstractDict{MOI.VariableIndex,MOI.VariableIndex} end
+Base.getindex(::IdentityMap, vi::MOI.VariableIndex) = vi
+
+function MOI.add_constraint(
+    model::MatrixOfConstraints{T},
+    func::MOI.AbstractFunction,
+    set::MOI.AbstractSet,
+) where {T}
+    i = set_index(model.sets, typeof(set))
+    if i === nothing || typeof(func) != _affine_function_type(T, typeof(set))
+        throw(MOI.UnsupportedConstraint{typeof(func),typeof(set)}())
+    end
+    return _add_constraint(model, i, IdentityMap(), func, set)
+end
+
+function _allocate_constraints(
+    model::MatrixOfConstraints{T},
+    src,
+    index_map,
+    ::Type{F},
+    ::Type{S},
+    filter_constraints::Union{Nothing,Function},
+) where {T,F,S}
+    i = set_index(model.sets, S)
+    if i === nothing || F != _affine_function_type(T, S)
+        throw(MOI.UnsupportedConstraint{F,S}())
+    end
+    cis_src = MOI.get(
+        src,
+        MOI.ListOfConstraintIndices{_affine_function_type(T, S),S}(),
+    )
+    if filter_constraints !== nothing
+        filter!(filter_constraints, cis_src)
+    end
+    for ci_src in cis_src
+        func = MOI.get(src, MOI.CanonicalConstraintFunction(), ci_src)
+        set = MOI.get(src, MOI.ConstraintSet(), ci_src)
+        push!(model.are_indices_mapped, length(model.caches) + 1)
+        index_map[ci_src] = _add_constraint(model, i, index_map, func, set)
+    end
+    return
+end
+
+function _load_constants(
+    constants,
+    offset,
+    func::MOI.AbstractScalarFunction,
+    set::MOI.AbstractScalarSet,
+)
+    MOI.throw_if_scalar_and_constant_not_zero(func, typeof(set))
+    return load_constants(constants, offset, set)
+end
+function _load_constants(
+    constants,
+    offset,
+    func::MOI.AbstractVectorFunction,
+    set::MOI.AbstractVectorSet,
+)
+    return load_constants(constants, offset, func)
+end
+
+function _load_constraints(
+    dest::MatrixOfConstraints,
+    index_map,
+    offset,
+    func_sets,
+)
+    for i in eachindex(func_sets)
+        func, set = func_sets[i]
+        if i in dest.are_indices_mapped
+            load_terms(dest.coefficients, index_map, func, offset)
+        else
+            load_terms(dest.coefficients, IdentityMap(), func, offset)
+        end
+        _load_constants(dest.constants, offset, func, set)
+        offset += MOI.output_dimension(func)
+    end
+    return offset
+end
+
+_add_variable(model::MatrixOfConstraints) = add_column(model.coefficients)
+
+function pass_nonvariable_constraints(
+    dest::MatrixOfConstraints,
+    src::MOI.ModelLike,
+    index_map::IndexMap,
+    constraint_types,
+    pass_cons = copy_constraints;
+    filter_constraints::Union{Nothing,Function} = nothing,
+)
+    for (F, S) in constraint_types
+        _allocate_constraints(dest, src, index_map, F, S, filter_constraints)
+    end
+    return
+end
+
+function final_touch(model::MatrixOfConstraints, index_map)
+    num_rows = MOI.dimension(model.sets)
+    resize!(model.constants, num_rows)
+    set_number_of_rows(model.coefficients, num_rows)
+
+    offset = 0
+    for cache in model.caches
+        offset = _load_constraints(model, index_map, offset, cache)
+    end
+
+    final_touch(model.coefficients)
+    empty!(model.are_indices_mapped)
+    empty!(model.caches)
+    return
+end
diff --git a/src/Utilities/mockoptimizer.jl b/src/Utilities/mockoptimizer.jl
index e4e35af863..817efd5d86 100644
--- a/src/Utilities/mockoptimizer.jl
+++ b/src/Utilities/mockoptimizer.jl
@@ -761,6 +761,9 @@ function MOI.supports_incremental_interface(
     return !mock.needs_allocate_load &&
            MOI.supports_incremental_interface(mock.inner_model, copy_names)
 end
+function final_touch(uf::MockOptimizer, index_map)
+    return final_touch(uf.inner_model, index_map)
+end
 
 # Allocate-Load Interface
 function supports_allocate_load(mock::MockOptimizer, copy_names::Bool)
diff --git a/src/Utilities/model.jl b/src/Utilities/model.jl
index 591a95ebfa..70ab05fd98 100644
--- a/src/Utilities/model.jl
+++ b/src/Utilities/model.jl
@@ -20,6 +20,16 @@ _no_lower_bound(::Type{T}) where {T<:AbstractFloat} = typemin(T)
 _no_upper_bound(::Type{T}) where {T} = zero(T)
 _no_upper_bound(::Type{T}) where {T<:AbstractFloat} = typemax(T)
 
+"""
+    function _add_variable end
+
+This is called by `AbstractModel` to inform the `constraints` field that a
+variable has been added. This is similar to
+[`MathOptInterface.add_variable`](@ref) except that it should return `nothing`.
+"""
+function _add_variable end
+
+function _add_variable(::Nothing) end
 function MOI.add_variable(model::AbstractModel{T}) where {T}
     vi = VI(model.num_variables_created += 1)
     push!(model.single_variable_mask, 0x0)
@@ -28,6 +38,7 @@ function MOI.add_variable(model::AbstractModel{T}) where {T}
     if model.variable_indices !== nothing
         push!(model.variable_indices, vi)
     end
+    _add_variable(model.constraints)
     return vi
 end
 
@@ -786,6 +797,9 @@ function MOI.copy_to(dest::AbstractModel, src::MOI.ModelLike; kws...)
     return automatic_copy_to(dest, src; kws...)
 end
 MOI.supports_incremental_interface(::AbstractModel, ::Bool) = true
+function final_touch(model::AbstractModel, index_map)
+    return final_touch(model.constraints, index_map)
+end
 
 # Allocate-Load Interface
 # Even if the model does not need it and use default_copy_to, it could be used
diff --git a/src/Utilities/product_of_sets.jl b/src/Utilities/product_of_sets.jl
new file mode 100644
index 0000000000..792552db09
--- /dev/null
+++ b/src/Utilities/product_of_sets.jl
@@ -0,0 +1,331 @@
+"""
+    abstract type ProductOfSets{T} end
+
+Represents a cartesian product of sets of given types.
+"""
+abstract type ProductOfSets{T} end
+
+"""
+    set_index(sets::ProductOfSets, ::Type{S}) where {S<:MOI.AbstractSet}
+
+Return an integer corresponding to the index of the set type in the list given
+by [`set_types`](@ref). If this set is not part of the list then it returns
+`nothing`.
+"""
+set_index(sets::ProductOfSets, ::Type{S}) where {S<:MOI.AbstractSet} = nothing
+
+"""
+    set_types(sets::ProductOfSets)
+
+Return the list of the types of the sets allowed in the cartesian product.
+"""
+function set_types end
+
+"""
+    add_set(sets::OrderedProductOfSets, i)
+
+Add a scalar set of type index `i`.
+
+    add_set(sets::OrderedProductOfSets, i, dim)
+
+Add a vector set of type index `i` and dimension `dim`.
+
+Both method return a unique id of the set that
+can be used to reference this set.
+"""
+function add_set end
+
+"""
+    indices(sets::OrderedProductOfSets, ci::MOI.ConstraintIndex)
+
+Return the indices in `1:MOI.dimension(sets)` corresponding to the set of id
+`ci.value`. For scalar sets, this return an integer and for vector sets, this
+return an `UnitRange`.
+"""
+function indices end
+
+function _sets_code(esc_name, T, type_def, set_types...)
+    code = Expr(:block, type_def)
+    esc_types = esc.(set_types)
+    for (i, esc_type) in enumerate(esc_types)
+        method = push!(
+            code.args,
+            :(
+                function $MOIU.set_index(
+                    ::$esc_name{$T},
+                    ::Type{$esc_type},
+                ) where {$T}
+                    return $i
+                end
+            ),
+        )
+    end
+    push!(code.args, :(function $MOIU.set_types(::$esc_name{$T}) where {$T}
+        return [$(esc_types...)]
+    end))
+    return code
+end
+
+"""
+    abstract type MixOfScalarSets{T} <: ProductOfSets{T} end
+
+Product of scalar sets in the order the constraints are added, mixing the
+constraints of different types.
+"""
+abstract type MixOfScalarSets{T} <: ProductOfSets{T} end
+macro mix_of_scalar_sets(name, set_types...)
+    esc_name = esc(name)
+    T = esc(:T)
+    type_def = :(struct $esc_name{$T} <: $MOIU.MixOfScalarSets{$T}
+        set_ids::Vector{Int}
+        function $esc_name{$T}() where {$T}
+            return new(Int[])
+        end
+    end)
+    return _sets_code(esc_name, T, type_def, set_types...)
+end
+
+MOI.is_empty(sets::MixOfScalarSets) = isempty(sets.set_ids)
+MOI.empty!(sets::MixOfScalarSets) = empty!(sets.set_ids)
+MOI.dimension(sets::MixOfScalarSets) = length(sets.set_ids)
+indices(::MixOfScalarSets, ci::MOI.ConstraintIndex) = ci.value
+function add_set(sets::MixOfScalarSets, i)
+    push!(sets.set_ids, i)
+    return length(sets.set_ids)
+end
+function MOI.get(
+    sets::MixOfScalarSets{T},
+    ::MOI.ListOfConstraintTypesPresent,
+) where {T}
+    present = Set(sets.set_ids)
+    return [
+        (_affine_function_type(T, S), S) for
+        S in set_types(sets) if set_index(sets, S) in present
+    ]
+end
+function MOI.get(
+    sets::MixOfScalarSets,
+    ::MOI.NumberOfConstraints{F,S},
+) where {F,S}
+    i = set_index(sets, S)
+    return count(isequal(i), sets.set_ids)
+end
+function MOI.get(
+    sets::MixOfScalarSets,
+    ::MOI.ListOfConstraintIndices{F,S},
+) where {F,S}
+    i = set_index(sets, S)
+    return LazyMap{MOI.ConstraintIndex{F,S}}(
+        j -> MOI.ConstraintIndex{F,S}(j),
+        Base.Iterators.Filter(
+            j -> sets.set_ids[j] == i,
+            eachindex(sets.set_ids),
+        ),
+    )
+end
+function MOI.is_valid(
+    sets::MixOfScalarSets,
+    ci::MOI.ConstraintIndex{F,S},
+) where {F,S}
+    i = set_index(sets, S)
+    return i !== nothing &&
+           ci.value in eachindex(sets.set_ids) &&
+           sets.set_ids[ci.value] == i
+end
+
+"""
+    abstract type OrderedProductOfSets{T} <: ProductOfSets{T} end
+
+Product of sets in the order the constraints are added, grouping the
+constraints of the same types contiguously.
+"""
+abstract type OrderedProductOfSets{T} <: ProductOfSets{T} end
+macro product_of_sets(name, set_types...)
+    esc_name = esc(name)
+    T = esc(:T)
+    type_def = :(
+        struct $esc_name{$T} <: $MOIU.OrderedProductOfSets{$T}
+            num_rows::Vector{Int}
+            dimension::Dict{Tuple{Int,Int},Int}
+            function $esc_name{$T}() where {$T}
+                return new(zeros(Int, $(1 + length(set_types))), Dict{Int,Int}())
+            end
+        end
+    )
+    return _sets_code(esc_name, T, type_def, set_types...)
+end
+
+"""
+    abstract type OrderedProductOfScalarSets{T} <: OrderedProductOfSets{T} end
+
+Same as [`OrderedProductOfSets`](@ref) except that all types are scalar sets,
+which allows a more efficient implementation.
+"""
+abstract type OrderedProductOfScalarSets{T} <: OrderedProductOfSets{T} end
+macro product_of_scalar_sets(name, set_types...)
+    esc_name = esc(name)
+    T = esc(:T)
+    type_def = :(struct $esc_name{$T} <: $MOIU.OrderedProductOfScalarSets{$T}
+        num_rows::Vector{Int}
+        function $esc_name{$T}() where {$T}
+            return new(zeros(Int, $(1 + length(set_types))))
+        end
+    end)
+    return _sets_code(esc_name, T, type_def, set_types...)
+end
+
+MOI.is_empty(sets::OrderedProductOfSets) = all(iszero, sets.num_rows)
+function MOI.empty!(sets::OrderedProductOfSets)
+    fill!(sets.num_rows, 0)
+    empty!(sets.dimension)
+    return
+end
+function MOI.empty!(sets::OrderedProductOfScalarSets)
+    fill!(sets.num_rows, 0)
+    return
+end
+
+function MOI.dimension(sets::OrderedProductOfSets)
+    for i in 3:length(sets.num_rows)
+        sets.num_rows[i] += sets.num_rows[i-1]
+    end
+    return sets.num_rows[end]
+end
+function indices(
+    sets::OrderedProductOfSets{T},
+    ci::MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S},
+) where {T,S}
+    i = set_index(sets, S)
+    return sets.num_rows[i] + ci.value + 1
+end
+function indices(
+    sets::OrderedProductOfSets{T},
+    ci::MOI.ConstraintIndex{MOI.VectorAffineFunction{T},S},
+) where {T,S}
+    i = set_index(sets, S)
+    return (sets.num_rows[i] + ci.value) .+ (1:sets.dimension[(i, ci.value)])
+end
+function add_set(sets::OrderedProductOfSets, i)
+    offset = sets.num_rows[i+1]
+    sets.num_rows[i+1] = offset + 1
+    return offset
+end
+function add_set(sets::OrderedProductOfSets, i, dim)
+    offset = sets.num_rows[i+1]
+    sets.num_rows[i+1] = offset + dim
+    sets.dimension[(i, offset)] = dim
+    return offset
+end
+function _num_indices(sets::OrderedProductOfSets, ::Type{S}) where {S}
+    i = set_index(sets, S)
+    return sets.num_rows[i+1] - sets.num_rows[i]
+end
+function MOI.get(
+    sets::OrderedProductOfSets{T},
+    ::MOI.ListOfConstraintTypesPresent,
+) where {T}
+    return [
+        (_affine_function_type(T, S), S) for
+        S in set_types(sets) if !iszero(_num_indices(sets, S))
+    ]
+end
+struct UnevenIterator
+    i::Int
+    start::Int
+    stop::Int
+    dimension::Dict{Tuple{Int,Int},Int}
+end
+Base.IteratorSize(::UnevenIterator) = Base.SizeUnknown()
+function Base.iterate(it::UnevenIterator, cur = it.start)
+    if cur >= it.stop
+        return nothing
+    else
+        return (cur, cur + it.dimension[(it.i, cur)])
+    end
+end
+function Base.in(x, it::UnevenIterator)
+    return x in it.start:(it.stop-1) && haskey(it.dimension, (it.i, x))
+end
+function _range_iterator(
+    ::OrderedProductOfSets{T},
+    ::Int,
+    start::Int,
+    stop::Int,
+    ::Type{MOI.ScalarAffineFunction{T}},
+) where {T}
+    return start:(stop-1)
+end
+function _range_iterator(
+    sets::OrderedProductOfSets{T},
+    i::Int,
+    start::Int,
+    stop::Int,
+    ::Type{MOI.VectorAffineFunction{T}},
+) where {T}
+    return UnevenIterator(i, start, stop, sets.dimension)
+end
+function _range(
+    sets::OrderedProductOfSets{T},
+    ::Type{F},
+    ::Type{S},
+) where {T,F,S}
+    i = set_index(sets, S)
+    if F != _affine_function_type(T, S) || i === nothing
+        return nothing
+    else
+        return _range_iterator(
+            sets,
+            i,
+            0,
+            sets.num_rows[i+1] - sets.num_rows[i],
+            F,
+        )
+    end
+end
+_length(r::UnitRange) = length(r)
+_length(r::UnevenIterator) = count(_ -> true, r)
+function MOI.get(
+    sets::OrderedProductOfSets,
+    ::MOI.NumberOfConstraints{F,S},
+) where {F,S}
+    r = _range(sets, F, S)
+    if r === nothing
+        return 0
+    else
+        return _length(r)
+    end
+end
+function _empty(
+    ::OrderedProductOfSets{T},
+    ::Type{<:MOI.ScalarAffineFunction},
+) where {T}
+    return 1:0
+end
+function _empty(
+    sets::OrderedProductOfSets{T},
+    ::Type{<:MOI.VectorAffineFunction},
+) where {T}
+    return UnevenIterator(1, 1, 0, sets.dimension)
+end
+
+function MOI.get(
+    sets::OrderedProductOfSets,
+    ::MOI.ListOfConstraintIndices{F,S},
+) where {F,S}
+    rows = _range(sets, F, S)
+    if rows === nothing
+        # Empty iterator
+        rows = _empty(sets, F)
+    end
+    return LazyMap{MOI.ConstraintIndex{F,S}}(
+        i -> MOI.ConstraintIndex{F,S}(i),
+        rows,
+    )
+end
+function MOI.is_valid(
+    sets::OrderedProductOfSets,
+    ci::MOI.ConstraintIndex{F,S},
+) where {F,S}
+    r = _range(sets, F, S)
+    return r !== nothing && ci.value in r
+end
diff --git a/src/Utilities/sparse_matrix.jl b/src/Utilities/sparse_matrix.jl
new file mode 100644
index 0000000000..35905fd40b
--- /dev/null
+++ b/src/Utilities/sparse_matrix.jl
@@ -0,0 +1,191 @@
+import SparseArrays
+
+"""
+    abstract type AbstractIndexing end
+
+Indexing to be used for storing the row and column indices of
+`MutableSparseMatrixCSC`. See [`ZeroBasedIndexing`](@ref) and
+[`OneBasedIndexing`](@ref).
+"""
+abstract type AbstractIndexing end
+
+"""
+    struct ZeroBasedIndexing <: AbstractIndexing end
+
+Zero-based indexing: the `i`th row or column has index `i - 1`. This is useful
+when the vectors of row and column indices need to be communicated to a
+library using zero-based indexing such as C libraries.
+"""
+struct ZeroBasedIndexing <: AbstractIndexing end
+
+"""
+    struct ZeroBasedIndexing <: AbstractIndexing end
+
+One-based indexing: the `i`th row or column has index `i`. This enables an
+allocation-free conversion of [`MutableSparseMatrixCSC`](@ref) to
+`SparseArrays.SparseMatrixCSC`.
+"""
+struct OneBasedIndexing <: AbstractIndexing end
+
+_first_index(::ZeroBasedIndexing) = 0
+_first_index(::OneBasedIndexing) = 1
+_shift(x, ::ZeroBasedIndexing, ::ZeroBasedIndexing) = x
+_shift(x::Integer, ::ZeroBasedIndexing, ::OneBasedIndexing) = x + 1
+_shift(x::Array{<:Integer}, ::ZeroBasedIndexing, ::OneBasedIndexing) = x .+ 1
+_shift(x::Integer, ::OneBasedIndexing, ::ZeroBasedIndexing) = x - 1
+_shift(x, ::OneBasedIndexing, ::OneBasedIndexing) = x
+
+"""
+    mutable struct MutableSparseMatrixCSC{Tv,Ti<:Integer,I<:AbstractIndexing}
+        indexing::I
+        m::Int
+        n::Int
+        colptr::Vector{Ti}
+        rowval::Vector{Ti}
+        nzval::Vector{Tv}
+    end
+
+Matrix type loading sparse matrices in the Compressed Sparse Column format.
+The indexing used is `indexing`, see [`AbstractIndexing`](@ref). The other
+fields have the same meaning than for `SparseArrays.SparseMatrixCSC` except
+that the indexing is different unless `indexing` is `OneBasedIndexing`.
+
+The matrix is loaded in 5 steps:
+1) `MOI.empty!` is called.
+2) `MOI.Utilities.add_column` and `MOI.Utilities.allocate_terms` are called in
+   any order.
+3) `MOI.Utilities.set_number_of_rows` is called.
+4) `MOI.Utilities.load_terms` is called for each affine function.
+5) `MOI.Utilities.final_touch` is called.
+"""
+mutable struct MutableSparseMatrixCSC{Tv,Ti<:Integer,I<:AbstractIndexing}
+    indexing::I
+    m::Int # Number of rows
+    n::Int # Number of columns
+    colptr::Vector{Ti}
+    rowval::Vector{Ti}
+    nzval::Vector{Tv}
+    function MutableSparseMatrixCSC{Tv,Ti,I}() where {Tv,Ti<:Integer,I}
+        return new{Tv,Ti,I}(I(), 0, 0, Ti[], Ti[], Tv[])
+    end
+end
+
+function MOI.empty!(A::MutableSparseMatrixCSC)
+    A.m = 0
+    A.n = 0
+    resize!(A.colptr, 1)
+    A.colptr[1] = 0
+    empty!(A.rowval)
+    return empty!(A.nzval)
+end
+
+"""
+    add_column(A::MutableSparseMatrixCSC)
+
+Add a column to the matrix `A`.
+"""
+function add_column(A::MutableSparseMatrixCSC)
+    A.n += 1
+    push!(A.colptr, 0)
+    return
+end
+
+"""
+    set_number_of_rows(coefficients, n)
+
+This function sets the number of rows to `coefficients`. This allows it
+to preallocate necessary datastructures before the data is loaded with
+[`load_terms`](@ref).
+"""
+function set_number_of_rows(A::MutableSparseMatrixCSC, num_rows)
+    A.m = num_rows
+    for i in 3:length(A.colptr)
+        A.colptr[i] += A.colptr[i-1]
+    end
+    resize!(A.rowval, A.colptr[end])
+    return resize!(A.nzval, A.colptr[end])
+end
+
+"""
+    final_touch(A::MutableSparseMatrixCSC)
+
+Informs the matrix `A` that all functions have been added with `load_terms`.
+No more modification is allowed unless `MOI.empty!` is called.
+"""
+function final_touch(A::MutableSparseMatrixCSC)
+    for i in length(A.colptr):-1:2
+        A.colptr[i] = _shift(A.colptr[i-1], ZeroBasedIndexing(), A.indexing)
+    end
+    return A.colptr[1] = _first_index(A.indexing)
+end
+
+_variable(term::MOI.ScalarAffineTerm) = term.variable
+_variable(term::MOI.VectorAffineTerm) = _variable(term.scalar_term)
+function _allocate_terms(colptr, index_map, terms)
+    for term in terms
+        colptr[index_map[_variable(term)].value+1] += 1
+    end
+end
+
+"""
+    allocate_terms(A::MutableSparseMatrixCSC, index_map, func)
+
+Informs `A` that the terms of the function `func` where the variable
+indices are mapped with `index_map` will be loaded with [`load_terms`](@ref).
+The function `func` should be canonicalized, see [`is_canonical`](@ref).
+"""
+function allocate_terms(A::MutableSparseMatrixCSC, index_map, func)
+    return _allocate_terms(A.colptr, index_map, func.terms)
+end
+
+_output_index(::MOI.ScalarAffineTerm) = 1
+_output_index(term::MOI.VectorAffineTerm) = term.output_index
+_coefficient(term::MOI.ScalarAffineTerm) = term.coefficient
+_coefficient(term::MOI.VectorAffineTerm) = _coefficient(term.scalar_term)
+function _load_terms(colptr, rowval, nzval, index_map, terms, offset)
+    for term in terms
+        ptr = colptr[index_map[_variable(term)].value] += 1
+        rowval[ptr] = offset + _output_index(term)
+        nzval[ptr] = _coefficient(term)
+    end
+end
+
+"""
+    load_terms(A::MutableSparseMatrixCSC, index_map, func, offset)
+
+Loads the terms of `func` to `A` mapping the variable indices with `index_map`.
+The `i`th dimension of `func` is loaded at the `(offset + i)`th row of `A`. The
+function should be allocated first with [`allocate_terms`](@ref). The function
+`func` should be canonicalized, see [`is_canonical`](@ref).
+"""
+function load_terms(A::MutableSparseMatrixCSC, index_map, func, offset)
+    return _load_terms(
+        A.colptr,
+        A.rowval,
+        A.nzval,
+        index_map,
+        func.terms,
+        _shift(offset, OneBasedIndexing(), A.indexing),
+    )
+end
+
+"""
+    Base.convert(::Type{SparseMatrixCSC{Tv, Ti}}, A::MutableSparseMatrixCSC{Tv, Ti, I}) where {Tv, Ti, I}
+
+Converts `A` to a `SparseMatrixCSC`. Note that the field `A.nzval` is **not
+copied** so if `A` is modified after the call of this function, it can still
+affect the value returned. Moreover, if `I` is `OneBasedIndexing`, `colptr`
+and `rowval` are not copied either, i.e., the conversion is allocation-free.
+"""
+function Base.convert(
+    ::Type{SparseArrays.SparseMatrixCSC{Tv,Ti}},
+    A::MutableSparseMatrixCSC{Tv,Ti},
+) where {Tv,Ti}
+    return SparseArrays.SparseMatrixCSC{Tv,Ti}(
+        A.m,
+        A.n,
+        _shift(A.colptr, A.indexing, OneBasedIndexing()),
+        _shift(A.rowval, A.indexing, OneBasedIndexing()),
+        A.nzval,
+    )
+end
diff --git a/src/Utilities/struct_of_constraints.jl b/src/Utilities/struct_of_constraints.jl
index ccbadd6ffc..619f551d0e 100644
--- a/src/Utilities/struct_of_constraints.jl
+++ b/src/Utilities/struct_of_constraints.jl
@@ -1,5 +1,9 @@
 abstract type StructOfConstraints <: MOI.ModelLike end
 
+function _add_variable(model::StructOfConstraints)
+    return broadcastcall(_add_variable, model)
+end
+
 function _throw_if_cannot_delete(model::StructOfConstraints, vis, fast_in_vis)
     broadcastcall(model) do constrs
         if constrs !== nothing
diff --git a/src/Utilities/universalfallback.jl b/src/Utilities/universalfallback.jl
index d12547f29c..d9bf3f95b7 100644
--- a/src/Utilities/universalfallback.jl
+++ b/src/Utilities/universalfallback.jl
@@ -132,6 +132,10 @@ function MOI.supports_incremental_interface(
     return MOI.supports_incremental_interface(uf.model, copy_names)
 end
 
+function final_touch(uf::UniversalFallback, index_map)
+    return final_touch(uf.model, index_map)
+end
+
 # References
 function MOI.is_valid(uf::UniversalFallback, idx::MOI.VariableIndex)
     return MOI.is_valid(uf.model, idx)
diff --git a/src/Utilities/vector_of_constraints.jl b/src/Utilities/vector_of_constraints.jl
index 35a9db648f..0373801163 100644
--- a/src/Utilities/vector_of_constraints.jl
+++ b/src/Utilities/vector_of_constraints.jl
@@ -145,6 +145,8 @@ function MOI.modify(
     return
 end
 
+function _add_variable(::VectorOfConstraints) end
+
 # Deletion of variables in vector of variables
 
 function _remove_variable(v::VectorOfConstraints, vi::MOI.VariableIndex)
diff --git a/test/Utilities/copy.jl b/test/Utilities/copy.jl
index 08883ee557..dc3755c0ba 100644
--- a/test/Utilities/copy.jl
+++ b/test/Utilities/copy.jl
@@ -653,6 +653,7 @@ struct OnlyCopyConstraints{F,S} <: MOI.ModelLike
         return new{F,S}(MOIU.VectorOfConstraints{F,S}())
     end
 end
+function MOI.Utilities._add_variable(::OnlyCopyConstraints) end
 MOI.empty!(model::OnlyCopyConstraints) = MOI.empty!(model.constraints)
 function MOI.supports_constraint(
     model::OnlyCopyConstraints,
diff --git a/test/Utilities/matrix_of_constraints.jl b/test/Utilities/matrix_of_constraints.jl
new file mode 100644
index 0000000000..33bffbbee6
--- /dev/null
+++ b/test/Utilities/matrix_of_constraints.jl
@@ -0,0 +1,335 @@
+using SparseArrays, Test
+import MathOptInterface
+const MOI = MathOptInterface
+const MOIT = MOI.Test
+const MOIU = MOI.Utilities
+
+function _test_matrix_equal(A::SparseMatrixCSC, B::SparseMatrixCSC)
+    @test A.m == B.m
+    @test A.n == B.n
+    @test A.nzval == B.nzval
+    @test A.rowval == B.rowval
+    @test A.colptr == B.colptr
+end
+function _test_matrix_equal(
+    A::MOIU.MutableSparseMatrixCSC{Tv,Ti,I},
+    B::SparseMatrixCSC,
+) where {Tv,Ti,I}
+    @test A.m == B.m
+    @test A.n == B.n
+    @test A.nzval == B.nzval
+    if I <: MOIU.OneBasedIndexing
+        @test A.rowval == B.rowval
+        @test A.colptr == B.colptr
+    else
+        @test A.rowval == B.rowval .- 1
+        @test A.colptr == B.colptr .- 1
+    end
+    sA = convert(typeof(B), A)
+    @test typeof(sA) == typeof(B)
+    return _test_matrix_equal(sA, B)
+end
+
+function _test(
+    query_test,
+    test,
+    ConstantsType::Type,
+    ProductOfSetsType::Type,
+    A,
+    b,
+    bridged::Bool,
+    Indexing,
+)
+    optimizer = MOIU.GenericOptimizer{
+        Float64,
+        MOIU.MatrixOfConstraints{
+            Float64,
+            MOIU.MutableSparseMatrixCSC{Float64,Int,Indexing},
+            ConstantsType,
+            ProductOfSetsType,
+        },
+    }()
+    _inner(model::MOI.Bridges.LazyBridgeOptimizer) = _inner(model.model)
+    _inner(model::MOI.Utilities.CachingOptimizer) = _inner(model.optimizer)
+    _inner(model::MOI.Utilities.MockOptimizer) = _inner(model.inner_model)
+    _inner(model::MOI.Utilities.UniversalFallback) = _inner(model.model)
+    _inner(model::MOI.Utilities.AbstractModel) = model
+    _A(model::MOIU.AbstractModel) = model.constraints.coefficients
+    _b(model::MOIU.AbstractModel) = model.constraints.constants
+    if bridged
+        optimizer = MOI.Bridges.full_bridge_optimizer(optimizer, Float64)
+    end
+    config = MOIT.TestConfig(solve = false, query_number_of_constraints = false)
+    test(optimizer, config)
+    MOI.Utilities.final_touch(optimizer, MOI.Utilities.IdentityMap())
+    _test_matrix_equal(_A(_inner(optimizer)), A)
+    @test _b(_inner(optimizer)) == b
+    query_test(_inner(optimizer))
+
+    MOI.empty!(optimizer)
+    model = MOIU.CachingOptimizer(MOIU.Model{Float64}(), optimizer)
+    model.mode = MOIU.MANUAL
+    MOIU.reset_optimizer(model)
+    @test MOIU.state(model) == MOIU.EMPTY_OPTIMIZER
+    test(model, config)
+    MOIU.attach_optimizer(model)
+    _test_matrix_equal(_A(_inner(model)), A)
+    @test _b(_inner(model)) == b
+    query_test(_inner(optimizer))
+
+    MOI.empty!(optimizer)
+    model = MOIU.CachingOptimizer(
+        MOIU.Model{Float64}(),
+        MOIU.MockOptimizer(MOIU.UniversalFallback(optimizer)),
+    )
+    model.mode = MOIU.MANUAL
+    MOIU.reset_optimizer(model)
+    @test MOIU.state(model) == MOIU.EMPTY_OPTIMIZER
+    test(model, config)
+    MOIU.attach_optimizer(model)
+    _test_matrix_equal(_A(_inner(model)), A)
+    @test _b(_inner(model)) == b
+    query_test(_inner(optimizer))
+
+    if !bridged
+        MOI.empty!(optimizer)
+        model = MOIU.CachingOptimizer(
+            MOIU.Model{Float64}(),
+            MOI.Bridges.full_bridge_optimizer(
+                MOIU.CachingOptimizer(optimizer, MOIU.MANUAL),
+                Float64,
+            ),
+        )
+        model.mode = MOIU.MANUAL
+        MOIU.reset_optimizer(model)
+        @test MOIU.state(model) == MOIU.EMPTY_OPTIMIZER
+        test(model, config)
+        MOIU.attach_optimizer(model)
+        inner = model.optimizer.model.model_cache
+        _test_matrix_equal(_A(inner), A)
+        @test _b(inner) == b
+        query_test(inner)
+
+        MOI.empty!(optimizer)
+        model = MOIU.CachingOptimizer(
+            MOIU.Model{Float64}(),
+            MOI.Bridges.full_bridge_optimizer(optimizer, Float64),
+        )
+        model.mode = MOIU.MANUAL
+        MOIU.reset_optimizer(model)
+        @test MOIU.state(model) == MOIU.EMPTY_OPTIMIZER
+        test(model, config)
+        MOIU.attach_optimizer(model)
+        _test_matrix_equal(_A(_inner(model)), A)
+        @test _b(_inner(model)) == b
+        query_test(_inner(optimizer))
+    end
+    return
+end
+
+function _lp(model, ::MOI.Test.TestConfig{T}) where {T}
+    MOI.empty!(model)
+    x = MOI.add_variable(model)
+    fx = one(T) * MOI.SingleVariable(x)
+    y = MOI.add_variable(model)
+    fy = one(T) * MOI.SingleVariable(y)
+    MOI.add_constraint(model, 3fx + 2fy, MOI.EqualTo(T(5)))
+    MOI.add_constraint(model, fx, MOI.GreaterThan(zero(T)))
+    MOI.add_constraint(model, -fy, MOI.LessThan(zero(T)))
+    return MOI.add_constraint(model, 5fx - 4fy, MOI.Interval(T(6), T(7)))
+end
+
+function matrix_lp(T::Type, ProductOfSetsType::Type)
+    optimizer = MOIU.GenericOptimizer{
+        T,
+        MOIU.MatrixOfConstraints{
+            T,
+            MOIU.MutableSparseMatrixCSC{T,Int,Indexing},
+            MOI.Utilities.Box{T},
+            ProductOfSetsType,
+        },
+    }()
+    return MOI.Utilities.final_touch(optimizer)
+end
+
+MOIU.@mix_of_scalar_sets(
+    MixLP,
+    MOI.EqualTo{T},
+    MOI.GreaterThan{T},
+    MOI.LessThan{T},
+    MOI.Interval{T},
+)
+MOIU.@product_of_scalar_sets(
+    OrdsLP,
+    MOI.EqualTo{T},
+    MOI.GreaterThan{T},
+    MOI.LessThan{T},
+    MOI.Interval{T},
+)
+MOIU.@product_of_sets(
+    OrdLP,
+    MOI.EqualTo{T},
+    MOI.GreaterThan{T},
+    MOI.LessThan{T},
+    MOI.Interval{T},
+)
+
+@testset "contlinear $Indexing" for Indexing in [
+    MOIU.OneBasedIndexing,
+    MOIU.ZeroBasedIndexing,
+]
+    A2 = sparse([1, 1], [1, 2], ones(2))
+    b2 = MOI.Utilities.Box([-Inf], [1.0])
+    Alp = sparse(
+        [1, 1, 2, 3, 4, 4],
+        [1, 2, 1, 2, 1, 2],
+        Float64[3, 2, 1, -1, 5, -4],
+    )
+    blp = MOI.Utilities.Box([5, 0, -Inf, 6], [5, Inf, 0, 7])
+    F = MOI.ScalarAffineFunction{Float64}
+    @testset "$SetType" for SetType in
+                            [MixLP{Float64}, OrdsLP{Float64}, OrdLP{Float64}]
+        _test(
+            MOIT.linear2test,
+            MOI.Utilities.Box{Float64},
+            SetType,
+            A2,
+            b2,
+            false,
+            Indexing,
+        ) do optimizer
+            S = MOI.LessThan{Float64}
+            @test [(F, S), (MOI.SingleVariable, MOI.GreaterThan{Float64})] ==
+                  MOI.get(optimizer, MOI.ListOfConstraintTypesPresent())
+            @test 1 == MOI.get(optimizer, MOI.NumberOfConstraints{F,S}())
+            cis = MOI.get(optimizer, MOI.ListOfConstraintIndices{F,S}())
+            @test 1 == length(collect(cis))
+            for ci in cis
+                @test MOI.is_valid(optimizer, ci)
+                @test 1 == MOI.Utilities.rows(optimizer.constraints, ci)
+            end
+        end
+        _test(
+            _lp,
+            MOI.Utilities.Box{Float64},
+            SetType,
+            Alp,
+            blp,
+            false,
+            Indexing,
+        ) do optimizer
+            con_types = [
+                (F, MOI.EqualTo{Float64}),
+                (F, MOI.GreaterThan{Float64}),
+                (F, MOI.LessThan{Float64}),
+                (F, MOI.Interval{Float64}),
+            ]
+            @test con_types ==
+                  MOI.get(optimizer, MOI.ListOfConstraintTypesPresent())
+            for i in eachindex(con_types)
+                F, S = con_types[i]
+                @test 1 == MOI.get(optimizer, MOI.NumberOfConstraints{F,S}())
+                cis = MOI.get(optimizer, MOI.ListOfConstraintIndices{F,S}())
+                @test 1 == length(collect(cis))
+                for ci in cis
+                    @test MOI.is_valid(optimizer, ci)
+                    @test i == MOI.Utilities.rows(optimizer.constraints, ci)
+                end
+            end
+        end
+    end
+end
+
+MOIU.@product_of_sets(Nonneg, MOI.Nonnegatives)
+MOIU.@product_of_sets(NonposNonneg, MOI.Nonpositives, MOI.Nonnegatives)
+MOIU.@product_of_sets(NonnegNonpos, MOI.Nonnegatives, MOI.Nonpositives)
+
+@testset "contconic $Indexing" for Indexing in [
+    MOIU.OneBasedIndexing,
+    MOIU.ZeroBasedIndexing,
+]
+    function _lin3_query(optimizer, con_types)
+        @test con_types ==
+              MOI.get(optimizer, MOI.ListOfConstraintTypesPresent())
+        k = 0
+        for (F, S) in con_types
+            function bad_types(F, S)
+                @test 0 == @inferred MOI.get(
+                    optimizer,
+                    MOI.NumberOfConstraints{F,S}(),
+                )
+                @test isempty(
+                    @inferred MOI.get(
+                        optimizer,
+                        MOI.ListOfConstraintIndices{F,S}(),
+                    )
+                )
+                @test !MOI.is_valid(optimizer, MOI.ConstraintIndex{F,S}(1))
+            end
+            BadF = MOI.ScalarAffineFunction{Int}
+            BadS = MOI.EqualTo{Int}
+            bad_types(F, BadS)
+            bad_types(BadF, S)
+            bad_types(BadF, BadS)
+            n = div(2, length(con_types))
+            @test n ==
+                  @inferred MOI.get(optimizer, MOI.NumberOfConstraints{F,S}())
+            cis = MOI.get(optimizer, MOI.ListOfConstraintIndices{F,S}())
+            @test n == length(collect(cis))
+            for ci in cis
+                @test MOI.is_valid(optimizer, ci)
+                @test !MOI.is_valid(optimizer, typeof(ci)(-1))
+                k += 1
+                @test k:k == MOI.Utilities.rows(optimizer.constraints, ci)
+            end
+        end
+    end
+    # We test here that the constraints are reordered by defining the sets
+    # in two different orders and check that it affects `b`.
+    A = sparse([1, 2], [1, 1], ones(2))
+    b = [-1.0, 1.0]
+    F = MOI.VectorAffineFunction{Float64}
+    _test(
+        MOIT.lin3test,
+        Vector{Float64},
+        NonnegNonpos{Float64},
+        A,
+        b,
+        false,
+        Indexing,
+    ) do optimizer
+        return _lin3_query(
+            optimizer,
+            [(F, MOI.Nonnegatives), (F, MOI.Nonpositives)],
+        )
+    end
+    b = [1.0, -1.0]
+    _test(
+        MOIT.lin3test,
+        Vector{Float64},
+        NonposNonneg{Float64},
+        A,
+        b,
+        false,
+        Indexing,
+    ) do optimizer
+        return _lin3_query(
+            optimizer,
+            [(F, MOI.Nonpositives), (F, MOI.Nonnegatives)],
+        )
+    end
+    # Here, we test that it works of some constraints are bridged but not all.
+    A = sparse([1, 2], [1, 1], [1.0, -1.0])
+    b = -ones(2)
+    _test(
+        MOIT.lin3test,
+        Vector{Float64},
+        Nonneg{Float64},
+        A,
+        b,
+        true,
+        Indexing,
+    ) do optimizer
+        return _lin3_query(optimizer, [(F, MOI.Nonnegatives)])
+    end
+end