Totally agree, I love using JuMP. At the moment, Iâ€™m just curious about model construction time, and there the Gurobi API seems quite fast. However, I noticed something odd about the problem constructions: they are giving/solving different problems?

For example, I modified your improved version of my code to return the solution (if one exists) and came across a case where the C API finds an optimal solution to (some) problem but `JuMP`

says the problem is infeasible, so I donâ€™t know what the C API is solvingâ€¦ Which makes me think thereâ€™s an issue in how Iâ€™m trying to use the C API?

```
julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
OS: Linux (x86_64-linux-gnu)
CPU: 16 Ă— Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, skylake-avx512)
Threads: 1 on 16 virtual cores
Environment:
LD_LIBRARY_PATH = :/home/myuser/software/gurobi1000/linux64/lib:/home/myuser/software/gurobi1000/linux64/lib
using Gurobi
using JuMP
using LinearAlgebra
using Random
function solve_lp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
start_time = time()
env_p = Ref{Ptr{Cvoid}}()
error = GRBloadenv(env_p, "lp.log")
env = env_p[]
error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 1)
@assert error == 0
model_p = Ref{Ptr{Cvoid}}()
error = GRBnewmodel(env, model_p, "lp", 0, C_NULL, C_NULL, C_NULL, C_NULL, C_NULL)
@assert error == 0
model = model_p[]
m, n = size(A)
error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)
@assert error == 0
error = GRBaddvars(
model, # model
n, # : numvars
0, # : numnz
C_NULL, # : *vbeg
C_NULL, # : *vind
C_NULL, # : *vval
c, # : *obj
zeros(Cdouble, n), # : *lb
fill(Inf, n), # : *ub
C_NULL, # : *vtype
C_NULL # : **varnames
)
@assert error == 0
for i in 1:m
error = GRBaddconstr(
model, # : *model
n, # : numnz
Cint.(0:(n-1)), # : *cind
A[i, :], # : *cval
# GRB_LESS_EQUAL, # : sense
GRB_EQUAL, # : sense
b[i], # : rhs
C_NULL, # : *constrname
)
@assert error == 0
end
println("model setup = $(time() - start_time)")
solve_time = @elapsed error = GRBoptimize(model)
println("model solve = $solve_time")
sol = ones(T, n)
ret = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
@assert ret == 0
objval = Ref{Cdouble}()
ret = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
@assert ret == 0
ret = GRBfreemodel(model)
@assert ret == 0
GRBfreeenv(env)
return sol, objval
end
function solve_lp_jump(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
start_time = time()
m, n = size(A)
M = Model(Gurobi.Optimizer)
set_silent(M)
@variable(M, x[1:n] .>= 0)
@objective(M, Max, c'*x)
@constraint(M, A*x .<= b)
println("model setup = $(time() - start_time)")
solve_time = @elapsed optimize!(M)
println("model solve = $(solve_time)")
return value.(x), JuMP.objective_value(M)
end
function solve_lp_jump2(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
start_time = time()
m, n = size(A)
model = direct_model(Gurobi.Optimizer())
set_silent(model)
set_string_names_on_creation(model, false)
@variable(model, x[1:n] >= 0)
@objective(model, Max, c' * x)
@constraint(model, A* x .<= b)
println("model setup = $(time() - start_time)")
solve_time = @elapsed optimize!(model)
println("model solve = $solve_time")
return value.(x), objective_value(model)
end
```

The problem:

```
julia> Random.seed!(1234)
TaskLocalRNG()
julia> m = 50
50
julia> n = 100
100
julia> A = randn(Cdouble, m, n);
julia> b = randn(Cdouble, m);
julia> c = randn(Cdouble, n);
julia> @time sol1, objval1 = solve_lp_direct(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.02570796012878418
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (linux64)
CPU model: Intel(R) Xeon(R) W-2145 CPU @ 3.70GHz, instruction set [SSE2|AVX|AVX2|AVX512]
Thread count: 8 physical cores, 16 logical processors, using up to 16 threads
Optimize a model with 50 rows, 100 columns and 5000 nonzeros
Model fingerprint: 0xd65f6397
Coefficient statistics:
Matrix range [1e-04, 5e+00]
Objective range [2e-02, 3e+00]
Bounds range [0e+00, 0e+00]
RHS range [5e-02, 2e+00]
Presolve time: 0.01s
Presolved: 50 rows, 100 columns, 5000 nonzeros
Iteration Objective Primal Inf. Dual Inf. Time
0 4.6543195e+31 8.349456e+31 4.654320e+01 0s
92 5.2160968e+00 0.000000e+00 0.000000e+00 0s
Solved in 92 iterations and 0.01 seconds (0.01 work units)
Optimal objective 5.216096848e+00
model solve = 0.030953253
0.069868 seconds (7.81 k allocations: 477.694 KiB, 18.47% compilation time)
julia> sum(sol1)
29.62203562836316
julia> @time sol2, objval2 = solve_lp_jump(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.0030510425567626953
model solve = 0.003450798
ERROR: Gurobi Error 10005: Unable to retrieve attribute 'ObjVal'
Stacktrace:
[1] _check_ret
@ ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:368 [inlined] [0/1940]
[2] get(model::Gurobi.Optimizer, attr::MathOptInterface.ObjectiveValue)
@ Gurobi ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:3097
[3] get(b::MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, attr::MathOptInterface.ObjectiveValue)
@ MathOptInterface.Bridges ~/.julia/packages/MathOptInterface/GIbN0/src/Bridges/bridge_optimizer.jl:1002
[4] _get_model_attribute(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, attr::MathOptInterface.ObjectiveValue)
@ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/GIbN0/src/Utilities/cachingoptimizer.jl:828
[5] get
@ ~/.julia/packages/MathOptInterface/GIbN0/src/Utilities/cachingoptimizer.jl:876 [inlined]
[6] _moi_get_result(model::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{Gurobi.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, args::MathOptInterface.ObjectiveValue)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:680
[7] get(model::Model, attr::MathOptInterface.ObjectiveValue)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:700
[8] objective_value(model::Model; result::Int64)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:54
[9] objective_value
@ ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:50 [inlined]
[10] solve_lp_jump(A::Matrix{Float64}, b::Vector{Float64}, c::Vector{Float64})
@ Main ./REPL[30]:12
[11] top-level scope
@ ./timing.jl:262 [inlined]
[12] top-level scope
@ ./REPL[55]:0
julia> @time sol3, objval3 = solve_lp_jump2(A, b, c);
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
model setup = 0.023720979690551758
model solve = 0.026847763
ERROR: Gurobi Error 10005: Unable to retrieve attribute 'ObjVal'
Stacktrace:
[1] _check_ret
@ ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:368 [inlined]
[2] get(model::Gurobi.Optimizer, attr::MathOptInterface.ObjectiveValue)
@ Gurobi ~/.julia/packages/Gurobi/FliRK/src/MOI_wrapper/MOI_wrapper.jl:3097
[3] _moi_get_result(model::Gurobi.Optimizer, args::MathOptInterface.ObjectiveValue)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:671
[4] get(model::Model, attr::MathOptInterface.ObjectiveValue)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/optimizer_interface.jl:700
[5] objective_value(model::Model; result::Int64)
@ JuMP ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:54
[6] objective_value
@ ~/.julia/packages/JuMP/Xw2ck/src/objective.jl:50 [inlined]
[7] solve_lp_jump2(A::Matrix{Float64}, b::Vector{Float64}, c::Vector{Float64})
@ Main ./REPL[31]:13
[8] top-level scope
@ ./timing.jl:262 [inlined]
[9] top-level scope
@ ./REPL[56]:0
julia>
```

Thanks! I will try comparing the above with the dense implementation from this example.