diff --git a/Project.toml b/Project.toml
index dcaf0f8..9738abb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,7 +14,7 @@ MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
 BinaryProvider = "~0.5"
 CEnum = "0.3"
 GLPK_jll = "~4.64.0"
-MathOptInterface = "~0.9.5"
+MathOptInterface = "~0.9.15"
 julia = "1"
 
 [extras]
diff --git a/perf/runbench.jl b/perf/runbench.jl
new file mode 100644
index 0000000..3ff6c32
--- /dev/null
+++ b/perf/runbench.jl
@@ -0,0 +1,156 @@
+using SparseArrays
+using LinearAlgebra
+using Random
+using GLPK
+using MathOptInterface
+# using ProfileView
+const MOI = MathOptInterface
+
+using TimerOutputs
+
+struct RandomLP
+    rows::Int
+    cols::Int
+    dens::Float64
+end
+
+function generate_moi_problem(model, At, b, c;
+    var_bounds = true, scalar = true)
+    cols, rows = size(At)
+    x = MOI.add_variables(model, cols)
+    A_cols = rowvals(At)
+    A_vals = nonzeros(At)
+    if var_bounds
+        for col in 1:cols
+            MOI.add_constraint(model, MOI.SingleVariable(x[col]),
+                MOI.LessThan(10.0))
+            MOI.add_constraint(model, MOI.SingleVariable(x[col]),
+                MOI.GreaterThan(-10.0))
+        end
+    end
+    if scalar
+        for row in 1:rows
+            MOI.add_constraint(model, MOI.ScalarAffineFunction(
+                [MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]]) for i in nzrange(At, row)], 0.0),
+                MOI.LessThan(b[row]))
+        end
+    else
+        for row in 1:rows
+            MOI.add_constraint(model, MOI.VectorAffineFunction(
+                [MOI.VectorAffineTerm(1, 
+                    MOI.ScalarAffineTerm(A_vals[i], x[A_cols[i]])
+                ) for i in nzrange(At, row)], [-b[row]]),
+                MOI.Nonpositives(1))
+        end
+    end
+    objective = MOI.ScalarAffineFunction(
+        [MOI.ScalarAffineTerm(c[i], x[i]) for i in findall(!iszero, c)],
+            0.0)
+    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), objective)
+    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
+    return x
+end
+
+function random_data(seed, data)
+
+    rows = data.rows
+    cols = data.cols
+    density = data.dens
+
+    p_neg_element = 0.0 # 0.25
+
+    rng = Random.MersenneTwister(seed)
+
+    f_A(r, n) = ifelse.(rand(r, n) .> p_neg_element, 1, -1) .* (15 .+ 30 .* rand(r, n))
+
+    At = sprand(rng, cols, rows, density, f_A)
+    b = 50 * rand(rng, rows)
+
+    # not using signs now
+    sign = ifelse.(rand(rng, rows) .> 0.2, 'L', 'G')
+
+    f_c(r, n) = 20 .* 2 .* (rand(r, n) .- 0.5)
+    c = sprand(rng, cols, 0.5, f_c)
+
+    return At, b, c
+end
+
+function bridged_cache_and_solver()
+    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+        MOI.Utilities.MANUAL), Float64)
+    GLPK_ = GLPK.Optimizer()
+    MOI.set(GLPK_, MOI.Silent(), true)
+    return model, GLPK_
+end
+function cache_and_solver()
+    model = MOI.Utilities.CachingOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+        MOI.Utilities.MANUAL)
+    GLPK_ = GLPK.Optimizer()
+    MOI.set(GLPK_, MOI.Silent(), true)
+    return model, GLPK_
+end
+function bridged_cached_solver()
+    model = MOI.Bridges.full_bridge_optimizer(MOI.Utilities.CachingOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+        GLPK.Optimizer()), Float64)
+    MOI.set(model, MOI.Silent(), true)
+    return model
+end
+function cached_solver()
+    model = MOI.Utilities.CachingOptimizer(
+        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
+        GLPK.Optimizer())
+    MOI.set(model, MOI.Silent(), true)
+    return model
+end
+
+function time_build_and_solve(to_build, to_solve, At, b, c, scalar = true)
+    @timeit "build" x = generate_moi_problem(to_build, At, b, c, scalar = scalar)
+    if to_build !== to_solve
+        @timeit "copy" MOI.copy_to(to_solve, to_build, copy_names = false)
+    end
+    MOI.set(to_solve, MOI.TimeLimitSec(), 0.0010)
+    @time @timeit "opt" MOI.optimize!(to_solve)
+    val = MOI.get(to_solve, MOI.SolveTime())
+    println(val)
+    @show MOI.get(to_solve, MOI.ObjectiveValue())
+    @show MOI.get(to_solve, MOI.TerminationStatus())
+end
+
+function solve_GLPK(seed, data; time_limit_sec=Inf)
+
+    reset_timer!()
+
+    @timeit "data"  At, b, c = random_data(1, data)
+    for i in 1:seed
+        # mod(i,5) == 0 && GC.gc()
+        GC.gc()
+        bridged_cache, pure_solver = bridged_cache_and_solver()
+        @timeit "bc + s" time_build_and_solve(bridged_cache, pure_solver, At, b, c)
+        
+        GC.gc()
+        cache, pure_solver2 = cache_and_solver()
+        @timeit "c + s" time_build_and_solve(cache, pure_solver2, At, b, c)
+        
+        GC.gc()
+        full_solver = bridged_cached_solver()
+        @timeit "bcs" time_build_and_solve(full_solver, full_solver, At, b, c)
+        
+        GC.gc()
+        full_solver = bridged_cached_solver()
+        @timeit "bcs + v" time_build_and_solve(full_solver, full_solver, At, b, c, false)
+        
+        GC.gc()
+        cache_solver = cached_solver()
+        @timeit "cs" time_build_and_solve(cache_solver, cache_solver, At, b, c)
+
+    end
+
+    print_timer()
+
+end
+
+solve_GLPK(2, RandomLP(11, 11, 0.5); time_limit_sec=5)
+solve_GLPK(20, RandomLP(10000, 10000, 0.005); time_limit_sec=5)
\ No newline at end of file
diff --git a/src/GLPK.jl b/src/GLPK.jl
index 32b9b4d..543164a 100644
--- a/src/GLPK.jl
+++ b/src/GLPK.jl
@@ -43,6 +43,7 @@ if !(v"4.64.0" <= _GLPK_VERSION <= v"4.64.0")
 end
 
 include("MOI_wrapper/MOI_wrapper.jl")
+include("MOI_wrapper/MOI_copy.jl")
 include("MOI_wrapper/MOI_callbacks.jl")
 include("MOI_wrapper/deprecated_constants.jl")
 
diff --git a/src/MOI_wrapper/MOI_copy.jl b/src/MOI_wrapper/MOI_copy.jl
new file mode 100644
index 0000000..4e246f2
--- /dev/null
+++ b/src/MOI_wrapper/MOI_copy.jl
@@ -0,0 +1,250 @@
+const DoubleDicts = MOI.Utilities.DoubleDicts
+
+function _add_bounds(::Vector{Float64}, ub, i, s::MOI.LessThan{Float64})
+    ub[i] = s.upper
+    return
+end
+
+function _add_bounds(lb, ::Vector{Float64}, i, s::MOI.GreaterThan{Float64})
+    lb[i] = s.lower
+    return
+end
+
+function _add_bounds(lb, ub, i, s::MOI.EqualTo{Float64})
+    lb[i], ub[i] = s.value, s.value
+    return
+end
+
+function _add_bounds(lb, ub, i, s::MOI.Interval{Float64})
+    lb[i], ub[i] = s.lower, s.upper
+    return
+end
+
+_bound_type(::Type{MOI.Interval{Float64}}) = INTERVAL
+_bound_type(::Type{MOI.EqualTo{Float64}}) = EQUAL_TO
+_bound_type(::Type{S}) where S = NONE
+
+function _extract_bound_data(
+    src, mapping, lb, ub, bound_type, s::Type{S}
+) where {S}
+    dict = DoubleDicts.with_type(mapping.conmap, MOI.SingleVariable, S)
+    type = _bound_type(s)
+    list = MOI.get(src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
+    add_sizehint!(dict, length(list))
+    for con_index in list
+        f = MOI.get(src, MOI.ConstraintFunction(), con_index)
+        s = MOI.get(src, MOI.ConstraintSet(), con_index)
+        column = mapping[f.variable].value
+        _add_bounds(lb, ub, column, s)
+        bound_type[column] = type
+        dict[con_index] = MOI.ConstraintIndex{MOI.SingleVariable, S}(column)
+    end
+    return
+end
+
+_add_type(type, i, ::MOI.Integer) = begin type[i] = INTEGER end
+_add_type(type, i, ::MOI.ZeroOne) = begin type[i] = BINARY end
+
+function _extract_type_data(src, mapping, var_type, ::Type{S}) where {S}
+    dict = DoubleDicts.with_type(mapping.conmap, MOI.SingleVariable, S)
+    list = MOI.get(src, MOI.ListOfConstraintIndices{MOI.SingleVariable, S}())
+    add_sizehint!(dict, length(list))
+    for con_index in list
+        f = MOI.get(src, MOI.ConstraintFunction(), con_index)
+        column = mapping[f.variable].value
+        _add_type(var_type, column, S())
+        dict[con_index] = MOI.ConstraintIndex{MOI.SingleVariable, S}(column)
+    end
+    return
+end
+
+function _init_index_map(src)
+    x_src = MOI.get(src, MOI.ListOfVariableIndices())
+    N = Cint(length(x_src))
+    is_contiguous = true
+    for x in x_src  # assuming all indexes are different
+        if !(1 <= x.value <= N)
+            is_contiguous = false
+        end
+    end
+    mapping = is_contiguous ? MOIU.IndexMap(N) : MOIU.IndexMap()
+    for (i, x) in enumerate(x_src)
+        mapping[x] = MOI.VariableIndex(i)
+    end
+    return N, mapping
+end
+
+_bounds2(s::MOI.GreaterThan{Float64}) = (s.lower, Inf)
+_bounds2(s::MOI.LessThan{Float64}) = (-Inf, s.upper)
+_bounds2(s::MOI.EqualTo{Float64}) = (s.value, s.value)
+_bounds2(s::MOI.Interval{Float64}) = (s.lower, s.upper)
+
+function add_sizehint!(vec, n)
+    len = length(vec)
+    return sizehint!(vec, len + n)
+end
+
+function _extract_row_data(src, mapping, lb, ub, I, J, V, ::Type{S}) where {S}
+    dict = mapping.conmap[MOI.ScalarAffineFunction{Float64}, S]
+    list = MOI.get(
+        src, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64}, S}()
+    )::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}}
+    N = length(list)
+    add_sizehint!(lb, N)
+    add_sizehint!(ub, N)
+    add_sizehint!(dict, N)
+    # first loop caches functions and counts terms to be added
+    n_terms = 0
+    function_cache = Array{MOI.ScalarAffineFunction{Float64}}(undef, N)
+    for i in 1:N
+        pre_function = MOI.get(src, MOI.ConstraintFunction(), list[i])
+        f = if MOIU.is_canonical(pre_function)
+            pre_function
+        else
+            # no duplicates are allowed in GLPK
+            MOIU.canonical(pre_function)
+        end
+        function_cache[i] = f
+        l, u = _bounds2(MOI.get(src, MOI.ConstraintSet(), list[i]))
+        push!(lb, l - f.constant)
+        push!(ub, u - f.constant)
+        n_terms += length(f.terms)
+    end
+
+    non_zeros = length(I)
+    row = non_zeros == 0 ? 1 : I[end] + 1
+    # resize + setindex is faster than sizehint! + push
+    # makes difference because I, J, V can be huge
+    resize!(I, non_zeros + n_terms)
+    resize!(J, non_zeros + n_terms)
+    resize!(V, non_zeros + n_terms)
+    for i in 1:N
+        f = function_cache[i]
+        for term in f.terms
+            non_zeros += 1
+            I[non_zeros] = row
+            J[non_zeros] = Cint(mapping[term.variable_index].value::Int64)
+            V[non_zeros] = term.coefficient
+        end
+        row += 1
+        ind = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, S}(row)
+        dict[list[i]] = ind
+    end
+    return
+end
+
+function test_data(src, dest)
+    for (F, S) in MOI.get(src, MOI.ListOfConstraints())
+        if !MOI.supports_constraint(dest, F, S)
+            throw(MOI.UnsupportedConstraint{F, S}("GLPK.Optimizer does not support constraints of type $F-in-$S."))
+        end
+    end
+    fobj_type = MOI.get(src, MOI.ObjectiveFunctionType())
+    if !MOI.supports(dest, MOI.ObjectiveFunction{fobj_type}())
+        throw(MOI.UnsupportedAttribute(MOI.ObjectiveFunction(fobj_type)))
+    end
+    return
+end
+
+function _add_all_variables(
+    model::Optimizer, N, lower, upper, bound_type, var_type
+)
+    glp_add_cols(model, N)
+    sizehint!(model.variable_info, N)
+    for i in 1:N
+        bound = get_moi_bound_type(lower[i], upper[i], bound_type)
+        # We started from empty model.variable_info, hence we assume ordering
+        index = CleverDicts.add_item(
+            model.variable_info,
+            VariableInfo(MOI.VariableIndex(i), i, bound, var_type[i]),
+        )
+        glp_bound_type = get_glp_bound_type(lower[i], upper[i])
+        glp_set_col_bnds(model, i, glp_bound_type, lower[i], upper[i])
+        if var_type[i] == BINARY
+            model.num_binaries += 1
+        end
+        if var_type[i] == INTEGER
+            model.num_integers += 1
+        end
+    end
+    return
+end
+
+function _add_all_constraints(dest::Optimizer, rl, ru, I, J, V)
+    n_constraints = length(rl)
+    glp_add_rows(dest, n_constraints)
+    glp_load_matrix(dest, length(I), offset(I), offset(J), offset(V))
+    sizehint!(dest.affine_constraint_info, n_constraints)
+    for i in 1:n_constraints
+        # assume ordered indexing
+        # assume no range constraints
+        if rl[i] == ru[i]
+            glp_set_row_bnds(dest, i, GLP_FX, rl[i], ru[i])
+            CleverDicts.add_item(
+                dest.affine_constraint_info,
+                ConstraintInfo(i, MOI.EqualTo{Float64}(rl[i])),
+            )
+        elseif ru[i] == Inf
+            glp_set_row_bnds(dest, i, GLP_LO, rl[i], GLP_DBL_MAX)
+            CleverDicts.add_item(
+                dest.affine_constraint_info,
+                ConstraintInfo(i, MOI.GreaterThan{Float64}(rl[i])),
+            )
+        else
+            glp_set_row_bnds(dest, i, GLP_UP, -GLP_DBL_MAX, ru[i])
+            CleverDicts.add_item(
+                dest.affine_constraint_info,
+                ConstraintInfo(i, MOI.LessThan{Float64}(ru[i])),
+            )
+        end
+    end
+    return
+end
+
+function MOI.copy_to(
+    dest::Optimizer, src::MOI.ModelLike; copy_names::Bool = false, kwargs...
+)
+    @assert MOI.is_empty(dest)
+    test_data(src, dest)
+    N, mapping = _init_index_map(src)
+    cl, cu = fill(-Inf, N), fill(Inf, N)
+    bound_type = fill(NONE, N)
+    var_type = fill(CONTINUOUS, N)
+    _extract_bound_data(src, mapping, cl, cu, bound_type, MOI.GreaterThan{Float64})
+    _extract_bound_data(src, mapping, cl, cu, bound_type, MOI.LessThan{Float64})
+    _extract_bound_data(src, mapping, cl, cu, bound_type, MOI.EqualTo{Float64})
+    _extract_bound_data(src, mapping, cl, cu, bound_type, MOI.Interval{Float64})
+    _extract_type_data(src, mapping, var_type, MOI.Integer)
+    _extract_type_data(src, mapping, var_type, MOI.ZeroOne)
+    _add_all_variables(dest, N, cl, cu, bound_type, var_type)
+    rl, ru, I, J, V = Float64[], Float64[], Cint[], Cint[], Float64[]
+    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.GreaterThan{Float64})
+    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.LessThan{Float64})
+    _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.EqualTo{Float64})
+    # range constraints not supported
+    # _extract_row_data(src, mapping, rl, ru, I, J, V, MOI.Interval{Float64})
+    _add_all_constraints(dest, rl, ru, I, J, V)
+    # Copy model attributes:
+    # obj function and sense are passed here
+    MOIU.pass_attributes(dest, src, copy_names, mapping)
+    variables = MOI.get(src, MOI.ListOfVariableIndices())
+    MOIU.pass_attributes(dest, src, copy_names, mapping, variables)
+    # TODO(odow): fix copy_names = false.
+    pass_constraint_attributes(dest, src, false, mapping)
+    return mapping
+end
+
+function pass_constraint_attributes(dest, src, copy_names, mapping)
+    ctr_types = MOI.get(src, MOI.ListOfConstraints())
+    for (F,S) in ctr_types
+        pass_constraint_attributes(dest, src, copy_names, mapping, F, S)
+    end
+    return
+end
+function pass_constraint_attributes(
+    dest, src, copy_names, mapping, ::Type{F}, ::Type{S}
+) where {F,S}
+    indices = MOI.get(src, MOI.ListOfConstraintIndices{F, S}())
+    MOIU.pass_attributes(dest, src, copy_names, mapping, indices)
+    return
+end
diff --git a/src/MOI_wrapper/MOI_wrapper.jl b/src/MOI_wrapper/MOI_wrapper.jl
index 4f401c2..67b9480 100644
--- a/src/MOI_wrapper/MOI_wrapper.jl
+++ b/src/MOI_wrapper/MOI_wrapper.jl
@@ -1,6 +1,7 @@
 import MathOptInterface
 
 const MOI = MathOptInterface
+const MOIU = MOI.Utilities
 const CleverDicts = MOI.Utilities.CleverDicts
 
 @enum(TypeEnum, CONTINUOUS, BINARY, INTEGER)
@@ -26,6 +27,10 @@ mutable struct VariableInfo
     function VariableInfo(index::MOI.VariableIndex, column::Int)
         return new(index, column, NONE, CONTINUOUS, "", "", "", "")
     end
+    function VariableInfo(index::MOI.VariableIndex, column::Int, bound::BoundEnum,
+        type::TypeEnum)
+        return new(index, column, bound, type, "", "", "", "")
+    end
 end
 
 struct ConstraintKey
@@ -39,6 +44,7 @@ mutable struct ConstraintInfo
     set::MOI.AbstractSet
     name::String
     ConstraintInfo(set) = new(0, set, "")
+    ConstraintInfo(row, set) = new(row, set, "")
 end
 
 # Dummy callback function for internal use only. Responsible for updating the
@@ -355,10 +361,6 @@ end
 
 MOI.Utilities.supports_default_copy_to(::Optimizer, ::Bool) = true
 
-function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kwargs...)
-    return MOI.Utilities.automatic_copy_to(dest, src; kwargs...)
-end
-
 function MOI.get(model::Optimizer, ::MOI.ListOfVariableAttributesSet)
     return MOI.AbstractVariableAttribute[MOI.VariableName()]
 end
@@ -758,6 +760,29 @@ end
 
 const GLP_DBL_MAX = prevfloat(Inf)
 
+function get_glp_bound_type(lower, upper)
+    if lower == upper
+        GLP_FX
+    elseif lower <= -GLP_DBL_MAX
+        upper >= GLP_DBL_MAX ? GLP_FR : GLP_UP
+    else
+        upper >= GLP_DBL_MAX ? GLP_LO : GLP_DB
+    end
+end
+function get_moi_bound_type(lower, upper, type)
+    if type == INTERVAL
+        INTERVAL
+    elseif type == EQUAL_TO
+        EQUAL_TO
+    elseif lower == upper
+        LESS_AND_GREATER_THAN
+    elseif lower <= -GLP_DBL_MAX
+        upper >= GLP_DBL_MAX ? NONE : LESS_THAN
+    else
+        upper >= GLP_DBL_MAX ? GREATER_THAN : LESS_AND_GREATER_THAN
+    end
+end
+
 function _set_variable_bound(
     model::Optimizer,
     column::Int,
@@ -770,18 +795,8 @@ function _set_variable_bound(
     if upper === nothing
         upper = glp_get_col_ub(model, column)
     end
-    bound_type = if lower == upper
-        GLP_FX
-    elseif lower <= -GLP_DBL_MAX
-        upper >= GLP_DBL_MAX ? GLP_FR : GLP_UP
-    else
-        upper >= GLP_DBL_MAX ? GLP_LO : GLP_DB
-    end
-    if upper < lower
-        glp_set_col_bnds(model, column, bound_type, lower, upper)
-    else
-        glp_set_col_bnds(model, column, bound_type, lower, upper)
-    end
+    bound_type = get_glp_bound_type(lower, upper)
+    glp_set_col_bnds(model, column, bound_type, lower, upper)
     return
 end
 
@@ -1563,7 +1578,7 @@ function MOI.get(model::Optimizer, attr::MOI.DualStatus)
     (status, _) = _get_status(model)
     if status == MOI.OPTIMAL
         return MOI.FEASIBLE_POINT
-    elseif status == MOI.INFEASIBLE || status == MOI.LOCALLY_INFEASIBLE
+    elseif status == MOI.INFEASIBLE
         if _certificates_potentially_available(model)
             return MOI.INFEASIBILITY_CERTIFICATE
         end
diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl
index 4028c73..3870e1b 100644
--- a/test/MOI_wrapper.jl
+++ b/test/MOI_wrapper.jl
@@ -42,7 +42,7 @@ const CONFIG = MOIT.TestConfig()
 end
 
 @testset "Linear tests" begin
-@testset "Default Solver"  begin
+    @testset "Default Solver"  begin
         MOIT.contlineartest(OPTIMIZER, MOIT.TestConfig(basis = true), [
             # This requires an infeasiblity certificate for a variable bound.
             "linear12",
@@ -462,3 +462,26 @@ end
     MOI.set(model, MOI.RawParameter("msg_lev"), GLPK.MSG_OFF)
     @test MOI.get(model, MOI.RawParameter("msg_lev")) == GLPK.GLP_MSG_OFF
 end
+
+@testset "MOI.copy" begin
+    dest = GLPK.Optimizer()
+    src = MOI.Utilities.Model{Float64}()
+    MOI.Utilities.loadfromstring!(src, """
+        variables: a, b, c, d
+        minobjective: a + b + c + d
+        c1: a >= 1.0
+        c2: b <= 2.0
+        c3: c == 3.0
+        c4: d in Interval(-4.0, 4.0)
+        c5: a in Integer()
+        c6: b in ZeroOne()
+        c7: a + b >= -1.1
+        c8: a + b <= 2.2
+        c8: c + d == 2.2
+    """)
+    MOI.copy_to(dest, src; copy_names = true)
+    v = MOI.get(dest, MOI.ListOfVariableIndices())
+    @test length(v) == 4
+    names = MOI.get.(dest, MOI.VariableName(), v)
+    @test names == ["a", "b", "c", "d"]
+end