# Slow model construction in JuMP

I find that when I am modeling a problem in JuMP, the model construction takes a long time. My first thought was that maybe I am writing inefficient code, and I was wondering:

1. if there might be any quick fixes to my JuMP code
2. if interfacing with MathOptInterface directly would be faster
3. if interfacing with Gurobi (or whatever solver I use) directly would be faster
using JuMP, Gurobi
m = 500
n = 100_000
A = randn(m, n)
b = randn(m)
lb = fill(-1.0, n)
ub = fill(+1.0, n)
lb[1] = -Inf
ub[2] = +Inf
julia> @time begin
model = Model(Gurobi.Optimizer)
@variable(model, x[1:n])
@variable(model, y[1:m])
@variable(model, z[1:m])
@variable(model, t)
@objective(model, Min, sum(x[i]^2 for i in 1:n)) # can replace with no objective and still slow
@constraint(model, y .== A * x .+ b)
@constraint(model, z .>= t .+ y)
@constraint(model, sum(z) .>= t)
@constraint(model, lb .<= x .<= ub)
end
GC: pause 67.38ms. collected 3342.926024MB. incr
GC: pause 396.65ms. collected 970.717568MB. full
GC: pause 429.67ms. collected 2750.147009MB. full
481.022531 seconds (3.94 M allocations: 12.238 GiB, 0.19% gc time, 0.02% compilation time)


There are a couple of minor changes you can do to speed things up (using n=10_000, not 100_000):

julia> GC.gc()

julia> @time begin
model = Model(Gurobi.Optimizer)
set_silent(model)
@variable(model, x[1:n])
@variable(model, y[1:m])
@variable(model, z[1:m])
@variable(model, t)
@objective(model, Min, sum(x[i]^2 for i in 1:n))
@constraint(model, y .== A * x .+ b)
@constraint(model, z .>= t .+ y)
@constraint(model, sum(z) .>= t)
@constraint(model, lb .<= x .<= ub)
optimize!(model)
end
7.960642 seconds (1.17 M allocations: 1.597 GiB, 1.20% gc time, 1.44% compilation time)

julia> GC.gc()

julia> @time begin
model = direct_model(Gurobi.Optimizer())
set_silent(model)
set_string_names_on_creation(model, false)
@variable(model, lb[i] <= x[i=1:n] <= ub[i])
@variable(model, y[1:m])
@variable(model, z[1:m])
@variable(model, t)
@objective(model, Min, sum(x[i]^2 for i in 1:n))
@constraint(model, [i=1:m], y[i] == sum(A[i, j] * x[j] for j in 1:n) + b[i])
@constraint(model, [i=1:m], z[i] >= t + y[i])
@constraint(model, sum(z) >= t)
optimize!(model)
end
7.030907 seconds (25.37 M allocations: 1.170 GiB, 2.17% gc time, 3.04% compilation time)


For example, @constraint(model, lb .<= x .<= ub) is actually a ScalarAffineFunction-in-Interval constraint, which JuMP needs to bridge. Youâ€™re much better off passing the bounds explicitly in @variable.

But youâ€™ll notice that it doesnâ€™t do very much.

The real problem is that your A matrix is dense.

500 * 100_000 is 50 million non-zeros. Gurobi is built with the assumption that the constraint matrix is sparse, typically with 1% or fewer non-zeros.

Something seems wrong with the formulation of your model though. It doesnâ€™t make sense. x = 0, and t = very large negative number is the optimal solution?

What is the problem you are actually trying to solve?

1 Like

Youâ€™re right this is not the real model. The problem Iâ€™m considering is a linear-reformulation of a CVaR constraint on y=Ax+b, and I was sloppy about the details here, though it has the main idea. As a first step, Iâ€™m testing some different solvers and was comparing them on randomly generated (dense) instances of different sizes.

Iâ€™m not at the point of testing on real problems yet with actual A so was using a dense A of a certain size as a rough (bad?) proxy for performance on even larger (but sparse) A matrices.

Do you think it will be faster to work directly in Gurobiâ€™s interface, i.e., GRBnewmodelâ€¦?

Using random unstructured dense data as a proxy isnâ€™t very useful. Real problems have sparse structure, and solvers are built to exploit that.

Going via the C API is probably a little faster, but youâ€™ll still have the same issue. You canâ€™t pass a dense matrix to Gurobi, you have to pass the row and column indices, as well as a vector of data. So you end up passing a a few 5e7 length vectors, and the data structures and algorithms inside Gurobi assume sparsity.

For most realistic problems, the cost of building the JuMP model is (much) less than the cost of solving it, so if youâ€™re struggling to build it, Gurobi will likely struggle to solve it.

You also have to factor in the time it takes you to write the C API version, as well as the likelihood of bugs etc . In general, you should probably just use JuMP, and only once all else fails should you consider alternatives.

2 Likes

I tried making/solving an LP from data A,b,c for the problem \max\{c^\top x : Ax \leq b, x \geq 0\} with the GRBnewmodel, GRBaddvars, and GRBaddconstr functions in the dense case (just for fun). To verify my construction, I compare against solving the same problem in JuMP. My code below shows that I am making an error in constructing the problem using the C interface. So I am asking Gurobi to see if they have an example of using the GRBnewmodel, GRBaddvars, and GRBaddconstr functions (or if they can diagnose my below code). Iâ€™m posting here as well just in case youâ€™re interested

For m=500 and n=10_000, it seems ~10x faster to use Gurobi directly.

#=
maximize    c' * x
subject to  A * x <= b
x >= 0
=#

using Gurobi
using JuMP

function solve_lp_direct(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
setup_time = @elapsed begin

# initialize model
start_time = time()
env_p = Ref{Ptr{Cvoid}}()
env = env_p[]
model_p = Ref{Ptr{Cvoid}}()
error = GRBnewmodel(
env,
model_p,
"lp",
0,
C_NULL,
C_NULL,
C_NULL,
C_NULL,
C_NULL
)
model = model_p[]
m, n = size(A)

# maximize
error = GRBsetintattr(model, GRB_INT_ATTR_MODELSENSE, GRB_MAXIMIZE)

# output
optimstatus = Ref{Cint}()
objval = Ref{Cdouble}()
sol = ones(Tf, n)

# variables and objective
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
)

# constraint
for i in 1:size(A,1)
Ai = A[i,:]
cind = findall(Ai .!= 0.0)
numnz = Cint(length(cind))
cval = Ai[cind]
cind .-= 1
cind = Cint.(cind)
sense = GRB_LESS_EQUAL
rhs = b[i]
model,  # : *model
numnz,  # : numnz
cind,   # : *cind
cval,   # : *cval
sense,  # : sense
rhs,    # : rhs
C_NULL, # : *constrname
)
@assert(error == 0)
end
end # begin
println("model setup = $(setup_time)") solve_time = @elapsed error = GRBoptimize(model) println("model solve =$(solve_time)")
# error = GRBwrite(model, "lp.lp");
error = GRBgetintattr(model, GRB_INT_ATTR_STATUS, optimstatus);
error = GRBgetdblattr(model, GRB_DBL_ATTR_OBJVAL, objval);
error = GRBgetdblattrarray(model, GRB_DBL_ATTR_X, 0, n, sol);
return sol
end

function solve_lp_jump(A::Matrix{Tf}, b::Vector{Tf}, c::Vector{Tf}) where Tf<:Cdouble
setup_time = @elapsed begin
m, n = size(A)
M = Model(Gurobi.Optimizer)
@variable(M, x[1:n] .>= 0)
@objective(M, Max, c'*x)
@constraint(M, A*x .<= b)
end
println("model setup = $setup_time") solve_time = optimize!(M) println("model solve =$(solve_time)")
return value.(x)
end

m = 50
n = 100
A = randn(Cdouble, m, n)
b = randn(Cdouble, m)
c = randn(Cdouble, n)
@time sol1 = solve_lp_direct(copy(A), copy(b), copy(c));
@time sol2 = solve_lp_jump(copy(A), copy(b), copy(c));
sum(abs.(sol1-sol2))


Let me take a look.

Itâ€™s probably because youâ€™re not freeing the Gurobi memory at the end. I could trigger the segfault with your original code, but not with my version below. Working with the C API from Julia can be difficult. Gurobi.jl does stuff to ensure that we donâ€™t free memory at the wrong time:

Hereâ€™s the code I ran:

using Gurobi
using JuMP

function solve_lp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
start_time = time()
env_p = Ref{Ptr{Cvoid}}()
env = env_p[]
error = GRBsetintparam(env, GRB_INT_PAR_OUTPUTFLAG, 0)
@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
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
model,              # : *model
n,                  # : numnz
Cint.(0:(n-1)),     # : *cind
A[i, :],            # : *cval
GRB_LESS_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")
ret = GRBfreemodel(model)
@assert ret == 0
GRBfreeenv(env)
return
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
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
end

m = 500
n = 10_000
A = randn(Cdouble, m, n);
b = randn(Cdouble, m);
c = randn(Cdouble, n);
@time sol1 = solve_lp_direct(A, b, c);
@time sol2 = solve_lp_jump(A, b, c);
@time sol3 = solve_lp_jump2(A, b, c);


and here are the results:

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.19083690643310547
model solve = 4.981149123
5.296278 seconds (489.36 k allocations: 86.229 MiB, 1.23% gc time, 2.05% compilation time)
model setup = 5.112233877182007
model solve = 10.832208031
16.768803 seconds (69.07 M allocations: 3.192 GiB, 6.36% gc time)
model setup = 2.414860963821411
model solve = 5.093770094
7.588198 seconds (43.45 M allocations: 1.522 GiB, 4.90% gc time, 1.67% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.1080930233001709
model solve = 4.50650111
4.683605 seconds (241.50 k allocations: 71.713 MiB, 1.21% compilation time)
model setup = 2.224179983139038
model solve = 4.767565328
7.054845 seconds (43.53 M allocations: 1.673 GiB, 4.95% gc time, 0.89% compilation time)
model setup = 2.6410491466522217
model solve = 4.550919043
7.258149 seconds (43.37 M allocations: 1.518 GiB, 9.25% gc time, 0.91% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.11534619331359863
model solve = 4.978835491
5.165873 seconds (241.50 k allocations: 71.713 MiB, 1.13% compilation time)
model setup = 2.328436851501465
model solve = 5.220872358
7.617613 seconds (43.51 M allocations: 1.673 GiB, 5.08% gc time, 0.89% compilation time)
model setup = 2.448287010192871
model solve = 5.000077219
7.517537 seconds (43.36 M allocations: 1.518 GiB, 5.13% gc time, 0.92% compilation time)

julia> include("/Users/oscar/Desktop/grb.jl")
model setup = 0.11165499687194824
model solve = 5.316515968
5.502262 seconds (241.50 k allocations: 71.713 MiB, 1.10% compilation time)
model setup = 2.561717987060547
model solve = 5.696499657
8.333657 seconds (43.54 M allocations: 1.673 GiB, 5.18% gc time, 0.90% compilation time)
model setup = 2.7369630336761475
model solve = 5.279443178
8.087859 seconds (43.38 M allocations: 1.518 GiB, 6.56% gc time, 0.88% compilation time)


If we ignore the first run because of compilation latency, then yeah, the model setup is much slower via JuMP. But that isnâ€™t surprising because the Gurobi.Optimizer object initializes a whole lot more than our C code. It prepares data structures at the JuMP level, creates a variety of dictionaries to map Cint indices to JuMP variables, normalizes the constraints, etc.

My bigger takeaway is that JuMP took 8 seconds. Gurobi C API took 5. So for the cost of 3 seconds we avoided all of the issues that youâ€™re encountering trying to get the C API to work.

p.s. I posted a link on your Gurobi portal post to this discussion. Note that Gurobi doesnâ€™t provide official support for Julia, so they might not be much help.

The C examples are here: https://www.gurobi.com/documentation/10.0/examples/c_examples.html

1 Like

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}}()
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
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
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)

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);
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);
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);
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.

            # GRB_LESS_EQUAL,     # : sense
GRB_EQUAL,          # : sense


They arenâ€™t the same problem. JuMP is <= and Gurobi is =.

Iâ€™m just curious about model construction time, and there the Gurobi API seems quite fast.

If you just want to measure construction time, set a small time limit of 0.01 seconds and call optimize!. It often isnâ€™t correct to measure the time building with actually calling optimize because solvers may lazily defer construction until the actual solve. (See [2206.03866] JuMP 1.0 for a discussion on how to benchmark.)

1 Like

Oh gosh I missed that as I was playing around with the code, wow just saw, sorry

1 Like

Hopefully the last thing. I would like to pass a sparse matrix to Gurobi. Gurobi uses CSR but SparseArrays is CSC. Iâ€™m sure that thereâ€™s a better way, but Iâ€™ve converted my sparse constraint matrix A to CSR with the SparseMatricesCSR package, but GRBXaddconstrs is complaining (gurobi error 1018 = â€śduplicated indicesâ€ť error). Is there anything obvious that Iâ€™m doing wrong?

using SparseMatricesCSR, Random, Gurobi
Random.seed!(1010101)
m = 1_000
n = 10_000
A = rand(Cdouble, m, n);
b = rand(Cdouble, m);
c = randn(Cdouble, n);
function solve_biglp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
Acsr = SparseMatricesCSR.SparseMatrixCSR(transpose(sparse(transpose(A))));
start_time = time()
env_p = Ref{Ptr{Cvoid}}()
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
model,             # model
Base.Csize_t(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
#         model,              # : *model
#         n,                  # : numnz
#         Cint.(0:(n-1)),     # : *cind
#         A[i, :],            # : *cval
#         GRB_LESS_EQUAL,     # : sense
#         b[i],               # : rhs
#         C_NULL,             # : *constrname
#     )
#     @assert error == 0
# end
consense = fill(GRB_LESS_EQUAL, Acsr.m)
conname = fill(C_NULL, Acsr.m)
conname = ["$i" for i in eachindex(b)] error = GRBXaddconstrs( model, # : *model Acsr.m, # : numconstrs Base.Csize_t(nnz(Acsr)), # : numnz Base.Csize_t.(Acsr.rowptr.-1), # : *cbeg Acsr.colval.-1, # : *cind Acsr.nzval, # : *cval consense, # : sense b, # : *rhs conname, # : **constrname ) @assert error == 0 GRBupdatemodel(model) GRBsetdblparam(env, "TimeLimit", 0.0) println("model setup =$(time() - start_time)")
solve_time = @elapsed error = GRBoptimize(model)
println("model setup = $(time() - start_time)") 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



When using the C API, you need to be very, very careful about the types. There are no error checks or warnings. Doing the wrong thing will just give an error, the wrong answer, or a segfault.

In this case, read the documentation:

The signature is

But your Acsr.colval.-1 is a Vector{Int}, when it needs to be Vector{Cint}. Gurobi will happily accept your argument, but it will reinterpret your Int as two packed Cints.

The fix is probably just Cint.(Acsr.colval.-1), although I havenâ€™t tested.

Your senses are probably wrong as well. I think it needs to be consense = fill(Cchar(GRB_LESS_EQUAL), Acsr.m).

1 Like

Thanks!! It worked with Cint and Cchar Just to document, the working code is below:

using Gurobi
using JuMP
using LinearAlgebra
using Random
using SparseMatricesCSR

function solve_sparselp_direct(A::Matrix{T}, b::Vector{T}, c::Vector{T}) where {T}
Acsr = SparseMatricesCSR.SparseMatrixCSR(transpose(sparse(transpose(A))));
start_time = time()
env_p = Ref{Ptr{Cvoid}}()
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
model,             # model
Base.Csize_t(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
consense = fill(Cchar(GRB_LESS_EQUAL), Acsr.m)
conname = fill(C_NULL, Acsr.m)
model,                         # : *model
Acsr.m,                        # : numconstrs
Base.Csize_t(nnz(Acsr)),       # : numnz
Base.Csize_t.(Acsr.rowptr.-1), # : *cbeg
Cint.(Acsr.colval.-1),         # : *cind
Acsr.nzval,                    # : *cval
consense,                      # : sense
b,                             # : *rhs
conname,                       # : **constrname
)
@assert error == 0
GRBupdatemodel(model)
GRBsetdblparam(env, "TimeLimit", 0.0)
println("model setup = $(time() - start_time)") solve_time = @elapsed error = GRBoptimize(model) println("model setup =$(time() - start_time)")
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

1 Like

Perhaps a clarification of what is actually going on will be helpful.

The definition of GRBXaddconstrs  is

function GRBXaddconstrs(model, numconstrs, numnz, cbeg, cind, cval, sense, rhs, constrnames)
ccall((:GRBXaddconstrs, libgurobi), Cint, (Ptr{GRBmodel}, Cint, Csize_t, Ptr{Csize_t}, Ptr{Cint}, Ptr{Cdouble}, Ptr{Cchar}, Ptr{Cdouble}, Ptr{Ptr{Cchar}}), model, numconstrs, numnz, cbeg, cind, cval, sense, rhs, constrnames)
end


Juliaâ€™s ccall converts arguments that are passed to it to the type given in the corresponding position in the tuple that is the third argument.

It converts arguments using Base.unsafe_convert(T, Base.cconvert(T, arg)) where T is the desired type and arg is the argument.

So take numconstrs. ccall expects a Cint, but you passed it Acsr.m which is a Int64.

julia> y = Int64(1)
1

julia> z = Base.unsafe_convert(Cint, Base.cconvert(Cint, y))
1

julia> typeof(z)
Int32


Bingo. That worked.

But *cind wanted Ptr{Cint}, and you passed it a Vector{Int64}.

julia> x = [1, 2, 3]
3-element Vector{Int64}:
1
2
3

julia> typeof(x)
Vector{Int64} (alias for Array{Int64, 1})

julia> Base.unsafe_convert(Ptr{Cint}, Base.cconvert(Ptr{Cint}, x))
Ptr{Int32} @0x0000000149a046e0

julia> pointer(x)
Ptr{Int64} @0x0000000149a046e0


In this case, Julia converted the pointer of x to the right type, but it didnâ€™t convert the underlying data to the element type Cint.

So as a general rule, if youâ€™re passing scalar, you can be fast and loose with the type and it should convert okay. But if youâ€™re passing a vector, you need to make sure the element type is the bit-exact type expected by the C library.

1 Like

Ahh, thanks! Really helpful to see the scalar / vector (pointer) distinction in action. I now see how this explains why it worked with â€śnumconstrsâ€ť set to Ascr.m::Int64 even though the C API wants Int32.

1 Like