diff --git a/Project.toml b/Project.toml
index fa7b579..5323233 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,11 +3,13 @@ uuid = "2f4eb8e6-3e35-4ae4-8c7a-f3d7d9bf20ed"
 version = "0.1.0"
 
 [deps]
+FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
+FillArrays = "1.15.0"
 MathOptInterface = "1"
 julia = "1.6"
 
diff --git a/src/MatrixOptInterface.jl b/src/MatrixOptInterface.jl
index 0fbc923..951fde2 100644
--- a/src/MatrixOptInterface.jl
+++ b/src/MatrixOptInterface.jl
@@ -6,6 +6,7 @@
 module MatrixOptInterface
 
 using LinearAlgebra, SparseArrays
+import FillArrays
 using MathOptInterface
 const MOI = MathOptInterface
 const MOIU = MathOptInterface.Utilities
@@ -67,6 +68,6 @@ end
 include("sparse_matrix.jl")
 include("conic_form.jl")
 include("matrix_input.jl")
-include("change_form.jl")
+#include("change_form.jl")
 
 end
diff --git a/src/change_form.jl b/src/change_form.jl
index 2b12b3d..c056437 100644
--- a/src/change_form.jl
+++ b/src/change_form.jl
@@ -15,21 +15,6 @@ function change_form(::Type{LPForm{T,AT,VT}}, lp::LPForm) where {T,AT,VT}
     )
 end
 
-function change_form(
-    ::Type{LPForm{T,AT,VT}},
-    lp::LPStandardForm{T},
-) where {T,AT,VT}
-    return LPForm{T,AT,VT}(
-        lp.sense,
-        lp.c,
-        lp.A,
-        lp.b,
-        lp.b,
-        fill(zero(T), length(lp.c)),
-        fill(typemax(T), length(lp.c)),
-    )
-end
-
 function change_form(
     ::Type{LPForm{T,AT,VT}},
     lp::LPGeometricForm{T},
@@ -125,31 +110,6 @@ function change_form(
     return change_form(LPGeometricForm{T,AT,VT}, temp_lp)
 end
 
-function change_form(
-    ::Type{LPStandardForm{T,AT,VT}},
-    lp::LPStandardForm,
-) where {T,AT,VT}
-    return LPStandardForm(lp.sense, lp.c, lp.A, lp.b)
-end
-
-function change_form(
-    ::Type{LPStandardForm{T,AT,VT}},
-    lp::LPGeometricForm{T},
-) where {T,AT,VT}
-    new_A = hcat(lp.A, -lp.A, AT(I, length(lp.b), length(lp.b)))
-    new_c = vcat(lp.c, -lp.c, fill(0.0, length(lp.b)))
-    return LPStandardForm{T,AT,VT}(lp.sense, new_c, new_A, copy(lp.b))
-end
-
-function change_form(
-    ::Type{LPStandardForm{T,AT,VT}},
-    lp::F,
-) where {T,AT,VT,F<:AbstractLPForm{T}}
-    temp_lp = change_form(LPForm{T,AT,VT}, lp)
-    new_lp = change_form(LPGeometricForm{T,AT,VT}, temp_lp)
-    return change_form(LPStandardForm{T,AT,VT}, new_lp)
-end
-
 function change_form(
     ::Type{LPSolverForm{T,AT,VT}},
     lp::LPSolverForm,
diff --git a/src/matrix_input.jl b/src/matrix_input.jl
index f0d5812..bf5209d 100644
--- a/src/matrix_input.jl
+++ b/src/matrix_input.jl
@@ -79,204 +79,234 @@ function MOI.get(
     return _dot(model.A[ci.value, :])
 end
 
-"""
-    LPStandardForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}
-
-Represents a problem of the form:
-```
-sense ⟨c, x⟩
-s.t.  A x == b
-      x ≥ 0
-```
-"""
-struct LPStandardForm{T,AT<:AbstractMatrix{T},VT<:AbstractVector{T}} <:
-       AbstractLPForm{T}
-    sense::MOI.OptimizationSense
-    c::VT
-    A::AT
-    b::VT
+# We cannot use MOI.Utilities.Hyperrectangle in case the vector type is not `Vector`
+struct Hyperrectangle{T,LT<:AbstractVector{T},UT<:AbstractVector{T}} <: MOI.Utilities.AbstractVectorBounds
+    lower::LT
+    upper::UT
 end
 
-function MOI.get(
-    ::LPStandardForm{T},
-    ::MOI.ListOfConstraintTypesPresent,
-) where {T}
-    return [
-        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
-        (MOI.VectorOfVariables, MOI.Nonnegatives),
-    ]
-end
+MOI.Utilities.function_constants(::Hyperrectangle{T}, row) where {T} = zero(T)
 
-const EQ{T} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}}
-
-function MOI.get(
-    model::LPStandardForm{T},
-    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.EqualTo{T}},
-) where {T}
-    # TODO return `collect` with MOI v0.9.15 (see https://github.com/jump-dev/MathOptInterface.jl/pull/1110)
-    return collect(MOIU.LazyMap{EQ{T}}(
-        i -> EQ{T}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `EQ{T}`
-        1:size(model.A, 1),
-    ))
-end
-
-function MOI.get(model::LPStandardForm, ::MOI.ConstraintSet, ci::EQ)
-    return MOI.EqualTo(model.b[ci.value])
-end
-
-const NONNEG = MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Nonnegatives}
-
-function MOI.get(
-    ::LPStandardForm,
-    ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,MOI.Nonnegatives},
-)
-    return [NONNEG(1)]
+# TODO specialize for SparseVector
+function linear_function(c::AbstractVector{T}) where {T}
+    return MOI.ScalarAffineFunction(
+        [MOI.ScalarAffineTerm(c[i], MOI.VariableIndex(i)) for i in eachindex(c)],
+        zero(T),
+    )
 end
 
-function MOI.get(model::LPStandardForm, ::MOI.ConstraintFunction, ci::NONNEG)
-    return MOI.VectorOfVariables(MOI.get(model, MOI.ListOfVariableIndices()))
+function linear_objective(sense::MOI.OptimizationSense, c::AbstractVector{T}) where {T}
+    model = MOI.Utilities.ObjectiveContainer{T}()
+    MOI.set(model, MOI.ObjectiveSense(), sense)
+    func = linear_function(c)
+    MOI.set(model, MOI.ObjectiveFunction{typeof(func)}(), func)
+    return model
 end
 
-function MOI.get(model::LPStandardForm, ::MOI.ConstraintSet, ci::NONNEG)
-    return MOI.Nonnegatives(MOI.get(model, MOI.NumberOfVariables()))
+# /!\ Type piracy
+function MOI.Utilities.extract_function(
+    A::Matrix{T},
+    row::Integer,
+    constant::T,
+) where {T}
+    func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], constant)
+    for col in axes(A, 2)
+        val = A[row, col]
+        if !iszero(val)
+            push!(
+                func.terms,
+                MOI.ScalarAffineTerm(val, MOI.VariableIndex(col)),
+            )
+        end
+    end
+    return func
 end
 
-"""
-    LPGeometricForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}} <: AbstractLPForm{T}
-
-Represents a linear problem of the form:
-```
-sense ⟨c, x⟩
-s.t.  Ax <= b
-```
-"""
-struct LPGeometricForm{T,AT<:AbstractMatrix{T},VT<:AbstractVector{T}} <:
-       AbstractLPForm{T}
-    sense::MOI.OptimizationSense
-    c::VT
-    A::AT
-    b::VT
+# Copy-paste of MOI.Utilities.MutableSparseMatrixCSC
+function _first_in_column(
+    A::SparseMatrixCSC,
+    row::Integer,
+    col::Integer,
+)
+    range = SparseArrays.nzrange(A, col)
+    idx = searchsortedfirst(view(A.rowval, range), row)
+    return get(range, idx, last(range) + 1)
 end
 
-function MOI.get(
-    ::LPGeometricForm{T},
-    ::MOI.ListOfConstraintTypesPresent,
+# Copy-paste of MOI.Utilities.MutableSparseMatrixCSC
+function MOI.Utilities.extract_function(
+    A::SparseMatrixCSC{T},
+    row::Integer,
+    constant::T,
 ) where {T}
-    return [(MOI.ScalarAffineFunction{T}, MOI.LessThan{T})]
+    func = MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{T}[], constant)
+    for col in 1:A.n
+        idx = _first_in_column(A, row, col)
+        if idx > last(SparseArrays.nzrange(A, col))
+            continue
+        end
+        r = A.rowval[idx]
+        if r == row
+            push!(
+                func.terms,
+                MOI.ScalarAffineTerm(A.nzval[idx], MOI.VariableIndex(col)),
+            )
+        end
+    end
+    return func
 end
 
-const LT{T} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},MOI.LessThan{T}}
-
-function MOI.get(
-    model::LPGeometricForm{T},
-    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},MOI.LessThan{T}},
-) where {T}
-    # FIXME `copy_constraint` needs a `Vector`
-    return collect(MOIU.LazyMap{LT{T}}(
-        i -> LT{T}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `LT{T}`
-        1:size(model.A, 1),
-    ))
+function free_variables(n, ::Type{T}) where {T}
+    model = MOI.Utilities.VariablesContainer{T}()
+    for _ in 1:n
+        MOI.add_variable(model)
+    end
+    return model
 end
 
-function MOI.get(model::LPGeometricForm, ::MOI.ConstraintSet, ci::LT)
-    return MOI.LessThan(model.b[ci.value])
+function nonnegative_variables(n, ::Type{T}) where {T}
+    model = MOI.Utilities.VariablesContainer{T}()
+    for _ in 1:n
+        vi = MOI.add_variable(model)
+        MOI.add_constraint(model, vi, MOI.GreaterThan(zero(T)))
+    end
+    return model
 end
 
-abstract type LPMixedForm{T} <: AbstractLPForm{T} end
-
-function MOI.get(
-    model::LPMixedForm{T},
-    ::MOI.ListOfConstraintTypesPresent,
-) where {T}
-    list = Tuple{DataType,DataType}[]
-    for S in
-        [MOI.EqualTo{T}, MOI.Interval{T}, MOI.GreaterThan{T}, MOI.LessThan{T}]
-        for F in [MOI.VariableIndex, MOI.ScalarAffineFunction{T}]
-            if !iszero(MOI.get(model, MOI.NumberOfConstraints{F,S}()))
-                push!(list, (F, S))
-            end
-        end
+function interval_variables(lower::AbstractVector{T}, upper::AbstractVector{T}) where {T}
+    @assert eachindex(lower) == eachindex(upper)
+    model = MOI.Utilities.VariablesContainer{T}()
+    for i in eachindex(lower)
+        vi = MOI.add_variable(model)
+        MOI.add_constraint(model, vi, MOI.GreaterThan(lower[i]))
+        MOI.add_constraint(model, vi, MOI.LessThan(upper[i]))
     end
-    return list
+    return model
 end
 
-struct BoundSense <: MOI.AbstractVariableAttribute end
+MOI.Utilities.@product_of_sets(
+    EqualTos,
+    MOI.EqualTo{T},
+)
 
-function MOI.get(model::LPMixedForm, ::BoundSense, vi::MOI.VariableIndex)
-    return _bound_sense(model.v_lb[vi.value], model.v_ub[vi.value])
+function equality_constraints(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
+    sets = EqualTos{T}()
+    for _ in eachindex(b)
+        MOI.Utilities.add_set(sets, MOI.Utilities.set_index(sets, MOI.EqualTo{T}))
+    end
+    MOI.Utilities.final_touch(sets)
+    constants = Hyperrectangle(b, b)
+    model = MOI.Utilities.MatrixOfConstraints{T}(A, constants, sets)
+    model.final_touch = true
+    return model
 end
 
-const LinearBounds{T} =
-    Union{MOI.EqualTo{T},MOI.Interval{T},MOI.GreaterThan{T},MOI.LessThan{T}}
+MOI.Utilities.@product_of_sets(
+    LessThans,
+    MOI.LessThan{T},
+)
 
-function MOI.get(
-    model::LPMixedForm{T},
-    ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T},S},
-) where {T,S<:LinearBounds{T}}
-    s = _sense(S)
-    return count(1:size(model.A, 1)) do i
-        return _constraint_bound_sense(model, i) == s
+function lessthan_constraints(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
+    sets = LessThans{T}()
+    for _ in eachindex(b)
+        MOI.Utilities.add_set(sets, MOI.Utilities.set_index(sets, MOI.LessThan{T}))
     end
+    MOI.Utilities.final_touch(sets)
+    constants = Hyperrectangle(FillArrays.Zeros{T}(length(b)), b)
+    model = MOI.Utilities.MatrixOfConstraints{T}(A, constants, sets)
+    model.final_touch = true
+    return model
 end
 
-const AFF{T,S} = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T},S}
-
-function MOI.get(
-    model::LPMixedForm{T},
-    ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T},S},
-) where {T,S<:LinearBounds{T}}
-    s = _sense(S)
-    # FIXME `copy_constraint` needs a `Vector`
-    return collect(
-        MOIU.LazyMap{AFF{T,S}}(
-            i -> AFF{T,S}(i), # FIXME `LazyMap` needs a `Function` so cannot just give `AFF{T, S}`
-            Base.Iterators.Filter(1:size(model.A, 1)) do i
-                return _constraint_bound_sense(model, i) == s
-            end,
-        ),
-    )
-end
+MOI.Utilities.@product_of_sets(
+    Intervals,
+    MOI.Interval{T},
+)
 
-const VBOUND{S} = MOI.ConstraintIndex{MOI.VariableIndex,S}
+function interval_constraints(A::AbstractMatrix{T}, lower::AbstractVector{T}, upper::AbstractVector{T}) where {T}
+    sets = Intervals{T}()
+    for _ in eachindex(lower)
+        MOI.Utilities.add_set(sets, MOI.Utilities.set_index(sets, MOI.Interval{T}))
+    end
+    MOI.Utilities.final_touch(sets)
+    constants = Hyperrectangle(lower, upper)
+    model = MOI.Utilities.MatrixOfConstraints{T}(A, constants, sets)
+    model.final_touch = true
+    return model
+end
+
+MOI.Utilities.@mix_of_scalar_sets(
+    MixedLinearSets,
+    MOI.EqualTo{T},
+    MOI.GreaterThan{T},
+    MOI.LessThan{T},
+)
 
-function MOI.get(
-    model::LPMixedForm{T},
-    ::MOI.NumberOfConstraints{MOI.VariableIndex,S},
-) where {T,S<:LinearBounds{T}}
-    s = _sense(S)
-    return count(MOI.get(model, MOI.ListOfVariableIndices())) do vi
-        return MOI.get(model, BoundSense(), vi) == s
+function mix_of_constraints(A::AbstractMatrix{T}, b::AbstractVector{T}, senses::Vector{ConstraintSense}) where {T}
+    @assert eachindex(b) == eachindex(senses)
+    sets = MixedLinearSets{T}()
+    for sense in senses
+        # Int(EQUAL_TO) is 0 but `MOI.Utilities.set_index(sets, MOI.EqualTo{T})` is 1
+        MOI.Utilities.add_set(sets, Int(sense) + 1)
     end
+    MOI.Utilities.final_touch(sets)
+    constants = Hyperrectangle(b, b)
+    model = MOI.Utilities.MatrixOfConstraints{T}(A, constants, sets)
+    model.final_touch = true
+    return model
 end
 
-function MOI.get(
-    model::LPMixedForm{T},
-    ::MOI.ListOfConstraintIndices{MOI.VariableIndex,S},
-) where {T,S<:LinearBounds{T}}
-    s = _sense(S)
-    return collect(
-        MOIU.LazyMap{VBOUND{S}}(
-            Base.Iterators.Filter(
-                MOI.get(model, MOI.ListOfVariableIndices()),
-            ) do vi
-                return MOI.get(model, BoundSense(), vi) == s
-            end,
-        ) do vi
-            return VBOUND{S}(vi.value)
-        end,
+"""
+    lp_standard_form(sense::MOI.OptimizationSense, c::AbstractVector, A::AbstractMatrix, b::AbstractVector)
+
+Represents a problem of the form:
+```
+sense ⟨c, x⟩
+s.t.  A x == b
+      x ≥ 0
+```
+"""
+function lp_standard_form(sense::MOI.OptimizationSense, c::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
+    m, n = size(A)
+    @assert length(c) == n
+    @assert length(b) == m
+    return MOI.Utilities.GenericModel{T}(
+        linear_objective(sense, c),
+        nonnegative_variables(n, T),
+        equality_constraints(A, b),
     )
 end
 
-function MOI.get(::LPMixedForm, ::MOI.ConstraintFunction, ci::VBOUND)
-    return MOI.VariableIndex(ci.value)
-end
+"""
+    lp_geometric_form(sense::MOI.OptimizationSense, c::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
 
-function MOI.get(model::LPMixedForm, ::MOI.ConstraintSet, ci::VBOUND)
-    return _bound_set(model.v_lb[ci.value], model.v_ub[ci.value])
+Represents a linear problem of the form:
+```
+sense ⟨c, x⟩
+s.t.  Ax <= b
+```
+"""
+function lp_geometric_form(sense::MOI.OptimizationSense, c::AbstractVector{T}, A::AbstractMatrix{T}, b::AbstractVector{T}) where {T}
+    m, n = size(A)
+    @assert length(c) == n
+    @assert length(b) == m
+    return MOI.Utilities.GenericModel{T}(
+        linear_objective(sense, c),
+        free_variables(n, T),
+        lessthan_constraints(A, b),
+    )
 end
 
 """
-    LPForm{T, AT<:AbstractMatrix{T}, VT <: AbstractVector{T}}
+    function lp_form(
+        sense::MOI.OptimizationSense,
+        c::AbstractVector{T},
+        A::AbstractMatrix{T},
+        c_lb::AbstractVector{T},
+        c_ub::AbstractVector{T},
+        v_lb::AbstractVector{T},
+        v_ub::AbstractVector{T},
+    ) where {T}
 
 Represents a problem of the form:
 ```
@@ -285,26 +315,38 @@ s.t.  c_lb <= Ax <= c_ub
       v_lb <=  x <= v_ub
 ```
 """
-struct LPForm{T,AT<:AbstractMatrix{T},VT<:AbstractVector{T}} <: LPMixedForm{T} #, V<:AbstractVector{T} #, M<:AbstractMatrix{T}}
-    sense::MOI.OptimizationSense
-    c::VT
-    A::AT
-    c_lb::VT
-    c_ub::VT
-    v_lb::VT
-    v_ub::VT
-end
-
-function _constraint_bound_sense(model::LPForm, i)
-    return _bound_sense(model.c_lb[i], model.c_ub[i])
-end
-
-function MOI.get(model::LPForm, ::MOI.ConstraintSet, ci::AFF)
-    return _bound_set(model.c_lb[ci.value], model.c_ub[ci.value])
+function lp_form(
+    sense::MOI.OptimizationSense,
+    c::AbstractVector{T},
+    A::AbstractMatrix{T},
+    c_lb::AbstractVector{T},
+    c_ub::AbstractVector{T},
+    v_lb::AbstractVector{T},
+    v_ub::AbstractVector{T},
+) where {T}
+    m, n = size(A)
+    @assert length(c) == n
+    @assert length(v_lb) == n
+    @assert length(v_ub) == n
+    @assert length(c_lb) == m
+    @assert length(c_ub) == m
+    return MOI.Utilities.GenericModel{T}(
+        linear_objective(sense, c),
+        interval_variables(v_lb, v_ub),
+        interval_constraints(A, c_lb, c_ub),
+    )
 end
 
 """
-    LPSolverForm{T, AT<:AbstractMatrix{T}, VT<:AbstractVector{T}}
+    lp_solver_form(
+        sense::MOI.OptimizationSense,
+        c::AbstractVector{T},
+        A::AbstractMatrix{T},
+        b::AbstractVector{T},
+        sense::Vector{ConstraintSense},
+        v_lb::AbstractVector{T},
+        v_ub::AbstractVector{T},
+    ) where {T}
 
 Represents a problem of the form:
 ```
@@ -313,32 +355,25 @@ s.t. Ax senses b
      v_lb <=  x <= v_ub
 ```
 """
-struct LPSolverForm{T,AT<:AbstractMatrix{T},VT<:AbstractVector{T}} <:
-       LPMixedForm{T}
-    sense::MOI.OptimizationSense
-    c::VT
-    A::AT
-    b::VT
-    senses::Vector{ConstraintSense}
-    v_lb::VT
-    v_ub::VT
-end
-
-function _constraint_bound_sense(model::LPSolverForm, i)
-    return model.senses[i]
-end
-
-function MOI.get(model::LPSolverForm, ::MOI.ConstraintSet, ci::AFF)
-    s = model.senses[ci.value]
-    β = model.b[ci.value]
-    if s == GREATER_THAN
-        return MOI.GreaterThan(β)
-    elseif s == LESS_THAN
-        return MOI.LessThan(β)
-    else
-        s == EQUAL_TO || error("Invalid $s")
-        return MOI.EqualTo(β)
-    end
+function lp_solver_form(
+    sense::MOI.OptimizationSense,
+    c::AbstractVector{T},
+    A::AbstractMatrix{T},
+    b::AbstractVector{T},
+    senses::Vector{ConstraintSense},
+    v_lb::AbstractVector{T},
+    v_ub::AbstractVector{T},
+) where {T}
+    m, n = size(A)
+    @assert length(c) == n
+    @assert length(v_lb) == n
+    @assert length(v_ub) == n
+    @assert length(b) == m
+    return MOI.Utilities.GenericModel{T}(
+        linear_objective(sense, c),
+        interval_variables(v_lb, v_ub),
+        mix_of_constraints(A, b, senses),
+    )
 end
 
 """
@@ -347,7 +382,7 @@ end
 A mixed-integer problem represented by a linear problem of type `LP`
 and a vector indicating each `VariableType`.
 """
-struct MILP{T,LP<:AbstractLPForm{T}}
+struct MILP{T,LP<:MOI.ModelLike}
     lp::LP
     variable_type::Vector{VariableType}
 end
diff --git a/test/linear.jl b/test/linear.jl
new file mode 100644
index 0000000..3eda2c1
--- /dev/null
+++ b/test/linear.jl
@@ -0,0 +1,189 @@
+module TestLinear
+
+using SparseArrays, Test
+import MathOptInterface as MOI
+import MatrixOptInterface as MatOI
+
+function runtests()
+    for name in names(@__MODULE__; all = true)
+        if startswith("$(name)", "test_")
+            @testset "$(name)" begin
+                getfield(@__MODULE__, name)()
+            end
+        end
+    end
+    return
+end
+
+function _test_expected(form, expected, var_names, con_names, S)
+    model = MOI.Utilities.Model{Float64}()
+    MOI.copy_to(
+        MOI.Bridges.Constraint.Scalarize{Float64}(model),
+        form,
+    )
+    MOI.set(
+        model,
+        MOI.VariableName(),
+        MOI.VariableIndex.(eachindex(var_names)),
+        var_names,
+    )
+    MOI.set(
+        model,
+        MOI.ConstraintName(),
+        MOI.ConstraintIndex{
+            MOI.ScalarAffineFunction{Float64},
+            S,
+        }.(eachindex(con_names)),
+        con_names,
+    )
+    return MOI.Test.util_test_models_equal(
+        model,
+        expected,
+        var_names,
+        con_names,
+    )
+end
+
+
+function test_standard_form()
+    dense_A = [
+        1.0 2.0
+        3.0 4.0
+    ]
+    dense_b = [5.0, 6.0]
+    dense_c = [7.0, 8.0]
+    s = """
+    variables: x1, x2
+    cx1: x1 >= 0.0
+    cx2: x2 >= 0.0
+    c1:  x1 + 2x2 == 5.0
+    c2: 3x1 + 4x2 == 6.0
+    minobjective: 7x1 + 8x2
+    """
+    expected = MOI.Utilities.Model{Float64}()
+    MOI.Utilities.loadfromstring!(expected, s)
+
+    var_names = ["x1", "x2"]
+    con_names = ["c1", "c2"]
+
+    sense = MOI.MIN_SENSE
+    v_lb = [0.0, 0.0]
+    v_ub = [Inf, Inf]
+
+    @testset "Matrix $(typeof(A))" for A in [dense_A, sparse(dense_A)]
+        @testset "Vector $(typeof(b))" for (b, c) in zip(
+            [dense_b, sparsevec(dense_b)],
+            [dense_c, sparsevec(dense_c)],
+        )
+            @testset "change $func" for (func, args) in [
+                (MatOI.lp_standard_form, (
+                    sense,
+                    c,
+                    A,
+                    b,
+                )),
+                (MatOI.lp_solver_form, (
+                    sense,
+                    c,
+                    A,
+                    b,
+                    [MatOI.EQUAL_TO, MatOI.EQUAL_TO],
+                    v_lb,
+                    v_ub,
+                )),
+            ]
+                lp = func(args...)
+                _test_expected(lp, expected, var_names, con_names, MOI.EqualTo{Float64})
+            end
+        end
+    end
+end
+
+function test_geometric_form()
+    dense_A = [
+        1.0 2.0
+        3.0 4.0
+    ]
+    dense_b = [5.0, 6.0]
+    dense_c = [7.0, 8.0]
+    s = """
+    variables: x1, x2
+    cx1: x1 >= 0.0
+    cx2: x2 >= 0.0
+    c1: x1 + 2x2 <= 5.0
+    c2: 3x1 + 4x2 <= 6.0
+    minobjective: 7x1 + 8x2
+    """
+    expected = MOI.Utilities.Model{Float64}()
+    MOI.Utilities.loadfromstring!(expected, s)
+
+    var_names = ["x1", "x2"]
+    con_names = ["c1", "c2"]
+
+    sense = MOI.MIN_SENSE
+
+    @testset "Matrix $(typeof(A))" for A in [dense_A, sparse(dense_A)]
+        @testset "Vector $(typeof(b))" for (b, c) in zip(
+            [dense_b, sparsevec(dense_b)],
+            [dense_c, sparsevec(dense_c)],
+        )
+            lp = MatOI.lp_geometric_form(
+                sense,
+                c,
+                A,
+                b,
+            )
+            _test_expected(lp, expected, var_names, con_names, MOI.LessThan{Float64})
+        end
+    end
+end
+
+function test_form()
+    dense_A = [
+        1.0 2.0
+        3.0 4.0
+    ]
+    dense_lb = [5.0, 6.0]
+    dense_ub = [7.0, 9.0]
+    dense_c = [7.0, 8.0]
+    s = """
+    variables: x1, x2
+    cx1: x1 >= 0.0
+    cx2: x2 >= 0.0
+    c1: x1 + 2x2 in MOI.Interval(5.0, 7.0)
+    c2: 3x1 + 4x2 in MOI.Interval(6.0, 9.0)
+    minobjective: 7x1 + 8x2
+    """
+    expected = MOI.Utilities.Model{Float64}()
+    MOI.Utilities.loadfromstring!(expected, s)
+
+    var_names = ["x1", "x2"]
+    con_names = ["c1", "c2"]
+
+    sense = MOI.MIN_SENSE
+    v_lb = [0.0, 0.0]
+    v_ub = [Inf, Inf]
+
+    @testset "Matrix $(typeof(A))" for A in [dense_A, sparse(dense_A)]
+        @testset "Vector $(typeof(c))" for (c_lb, c_ub, c) in zip(
+            [dense_lb, sparsevec(dense_lb)],
+            [dense_ub, sparsevec(dense_ub)],
+            [dense_c, sparsevec(dense_c)],
+        )
+            lp = MatOI.lp_form(
+                sense,
+                c,
+                A,
+                c_lb,
+                c_ub,
+                v_lb,
+                v_ub,
+            )
+            _test_expected(lp, expected, var_names, con_names, MOI.Interval{Float64})
+        end
+    end
+end
+
+end
+
+TestLinear.runtests()
diff --git a/test/runtests.jl b/test/runtests.jl
index 308ff0a..d95aa47 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -16,182 +16,5 @@ const MOIT = MOI.Test
 const ATOL = 1e-4
 const RTOL = 1e-4
 
+include("linear.jl")
 include("conic_form.jl")
-
-const dense_A = [
-    1.0 2.0
-    3.0 4.0
-]
-@testset "Matrix $(typeof(A))" for A in [dense_A, sparse(dense_A)]
-    dense_b = [5.0, 6.0]
-    dense_c = [7.0, 8.0]
-    @testset "Vector $(typeof(b))" for (b, c) in zip(
-        [dense_b, sparsevec(dense_b)],
-        [dense_c, sparsevec(dense_c)],
-    )
-        @testset "Standard form LP" begin
-            s = """
-            variables: x1, x2
-            cx1: x1 >= 0.0
-            cx2: x2 >= 0.0
-            c1:  x1 + 2x2 == 5.0
-            c2: 3x1 + 4x2 == 6.0
-            minobjective: 7x1 + 8x2
-            """
-            expected = MOIU.Model{Float64}()
-            MOIU.loadfromstring!(expected, s)
-
-            var_names = ["x1", "x2"]
-            con_names = ["c1", "c2"]
-            vcon_names = ["cx1", "cx2"]
-            model = MOIU.Model{Float64}()
-
-            sense = MOI.MIN_SENSE
-            v_lb = [0.0, 0.0]
-            v_ub = [Inf, Inf]
-
-            function test_expected(form)
-                MOI.copy_to(
-                    MOI.Bridges.Constraint.Scalarize{Float64}(model),
-                    form,
-                )
-                MOI.set(
-                    model,
-                    MOI.VariableName(),
-                    MOI.VariableIndex.(1:2),
-                    var_names,
-                )
-                MOI.set(
-                    model,
-                    MOI.ConstraintName(),
-                    MOI.ConstraintIndex{
-                        MOI.ScalarAffineFunction{Float64},
-                        MOI.EqualTo{Float64},
-                    }.(1:2),
-                    con_names,
-                )
-                return MOIT.util_test_models_equal(
-                    model,
-                    expected,
-                    var_names,
-                    con_names,
-                )
-            end
-
-            @testset "change $(typeof(lp))" for lp in [
-                MatOI.LPStandardForm{Float64,typeof(A),typeof(c)}(
-                    sense,
-                    c,
-                    A,
-                    b,
-                ),
-                MatOI.LPForm{Float64,typeof(A),typeof(c)}(
-                    sense,
-                    c,
-                    A,
-                    b,
-                    b,
-                    v_lb,
-                    v_ub,
-                ),
-                MatOI.LPSolverForm{Float64,typeof(A),typeof(c)}(
-                    sense,
-                    c,
-                    A,
-                    b,
-                    [MatOI.EQUAL_TO, MatOI.EQUAL_TO],
-                    v_lb,
-                    v_ub,
-                ),
-            ]
-                test_expected(lp)
-                @testset "to $F" for F in [
-                    #MatOI.LPStandardForm{Float64, typeof(A)}, # FIXME doesn't work as the form gets bloated in the conversion
-                    MatOI.LPForm{Float64,typeof(A),typeof(c)},
-                    MatOI.LPSolverForm{Float64,typeof(A),typeof(c)},
-                ]
-                    test_expected(MatOI.change_form(F, lp))
-                end
-            end
-        end
-        @testset "Geometric form LP" begin
-            s = """
-            variables: x1, x2
-            c1:  x1 + 3x2 <= 7.0
-            c2: 2x1 + 4x2 <= 8.0
-            maxobjective: 5x1 + 6x2
-            """
-            expected = MOIU.Model{Float64}()
-            MOIU.loadfromstring!(expected, s)
-
-            var_names = ["x1", "x2"]
-            con_names = ["c1", "c2"]
-            model = MOIU.Model{Float64}()
-
-            sense = MOI.MAX_SENSE
-            v_lb = [-Inf, -Inf]
-            v_ub = [Inf, Inf]
-
-            function test_expected(form)
-                MOI.copy_to(model, form)
-                MOI.set(
-                    model,
-                    MOI.VariableName(),
-                    MOI.VariableIndex.(1:2),
-                    var_names,
-                )
-                MOI.set(
-                    model,
-                    MOI.ConstraintName(),
-                    MOI.ConstraintIndex{
-                        MOI.ScalarAffineFunction{Float64},
-                        MOI.LessThan{Float64},
-                    }.(1:2),
-                    con_names,
-                )
-                return MOIT.util_test_models_equal(
-                    model,
-                    expected,
-                    var_names,
-                    con_names,
-                )
-            end
-
-            @testset "change $(typeof(lp))" for lp in [
-                MatOI.LPGeometricForm{Float64,typeof(A'),typeof(c)}(
-                    sense,
-                    b,
-                    A',
-                    c,
-                ),
-                MatOI.LPForm{Float64,typeof(A'),typeof(c)}(
-                    sense,
-                    b,
-                    A',
-                    [-Inf, -Inf],
-                    c,
-                    v_lb,
-                    v_ub,
-                ),
-                MatOI.LPSolverForm{Float64,typeof(A'),typeof(c)}(
-                    sense,
-                    b,
-                    A',
-                    c,
-                    [MatOI.LESS_THAN, MatOI.LESS_THAN],
-                    v_lb,
-                    v_ub,
-                ),
-            ]
-                test_expected(lp)
-                @testset "to $F" for F in [
-                    MatOI.LPGeometricForm{Float64,typeof(A),typeof(c)},
-                    MatOI.LPForm{Float64,typeof(A),typeof(c)},
-                    MatOI.LPSolverForm{Float64,typeof(A),typeof(c)},
-                ]
-                    test_expected(MatOI.change_form(F, lp))
-                end
-            end
-        end
-    end
-end