# Slow model construction in JuMP

**URL:** <https://discourse.julialang.org/t/slow-model-construction-in-jump/95455>\
**Category:** Optimization (Mathematical)\
**Tags:** jump\
**Created:** [March 2, 2023, 4:45pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455 "2023-03-02T16:45:54Z")\
**Posts on this page:** 15\
**Page:** 1

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 2, 2023, 4:45pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/1 "2023-03-02T16:45:54Z")

</div>

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

```julia
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                                                                                                                                          
Set parameter Username
Academic license - for non-commercial use only - expires 2023-12-17
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)

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 2, 2023, 8:22pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/2 "2023-03-02T20:22:42Z")

</div>

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

```plaintext
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?

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 2, 2023, 10:45pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/3 "2023-03-02T22:45:11Z")

</div>

Thanks again for your comments:)

> [@odow](#):
>
> 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?

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.

> [@odow](#):
>
> `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.

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`…?

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 2, 2023, 11:06pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/4 "2023-03-02T23:06:01Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 3, 2023, 1:30am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/5 "2023-03-03T01:30:50Z")

</div>

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](https://support.gurobi.com/hc/en-us/community/posts/13424957892241-Construct-LP-from-matrix-using-GRBnewmodel-GRBaddvars-and-GRBaddconstr) 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.

```julia
#=
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}}()
    error = GRBloadenv(env_p, "lp.log")
    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
    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
    )

    # 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]
      error = GRBaddconstr(
        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))

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 3, 2023, 1:33am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/6 "2023-03-03T01:33:26Z")

</div>

Let me take a look.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 3, 2023, 2:10am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/7 "2023-03-03T02:10:28Z")

</div>

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:

> <https://github.com/jump-dev/Gurobi.jl/blob/8c98d0d9988c67005bddfc2cf5b0365357551e5d/src/MOI_wrapper/MOI_wrapper.jl#L120-L131>

Here’s the code I ran:

```plaintext
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}}()
    error = GRBloadenv(env_p, "lp.log")
    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
    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
            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
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](https://www.gurobi.com/documentation/10.0/examples/c_examples.html)

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 3, 2023, 5:39am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/8 "2023-03-03T05:39:35Z")

</div>

> [@odow](#):
>
> 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.

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
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
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>

```

> [@odow](#):
>
> The C examples are here: [https://www.gurobi.com/documentation/10.0/examples/c\_examples.html](https://www.gurobi.com/documentation/10.0/examples/c_examples.html)

Thanks! I will try comparing the above with the dense implementation from [this](https://www.gurobi.com/documentation/10.0/examples/dense_c_c.html) example.

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 3, 2023, 5:52am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/9 "2023-03-03T05:52:23Z")

</div>

```julia
            # 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: Recent improvements to a modeling language for mathematical optimization](https://arxiv.org/abs/2206.03866) for a discussion on how to benchmark.)

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 3, 2023, 5:57am UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/10 "2023-03-03T05:57:47Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 14, 2023, 8:25pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/11 "2023-03-14T20:25:50Z")

</div>

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?

```julia
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}}()
    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 = GRBXaddvars(
        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
    # error = GRBaddconstr(
    # 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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 14, 2023, 9:21pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/12 "2023-03-14T21:21:06Z")

</div>

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:

- [https://www.gurobi.com/documentation/9.5/refman/c\_xaddconstrs.html](https://www.gurobi.com/documentation/9.5/refman/c_xaddconstrs.html)

The signature is

![image](https://global.discourse-cdn.com/julialang/original/3X/b/7/b7853b22cb2593408d5b81a2f982570e2a1dee86.png)

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 `Cint`s.

The fix is probably just `Cint.(Acsr.colval.-1)`, although I haven’t tested.

You should also follow this advice:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/7/0/704a0d74699fc2cef9e0d33626af4f580468c1bb.png)

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

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 14, 2023, 9:59pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/13 "2023-03-14T21:59:45Z")

</div>

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

```julia
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}}()
    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 = GRBXaddvars(
        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)
    error = GRBXaddconstrs(
        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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [March 14, 2023, 10:13pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/14 "2023-03-14T22:13:49Z")

</div>

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

The definition of `GRBXaddconstrs ` is

```julia
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
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
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.

---

<div class="post-metadata">

**Author:** ![jacob-roth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jacob-roth/32/1862_2.png) [@jacob-roth](https://discourse.julialang.org/u/jacob-roth)\
**Post date:** [March 14, 2023, 10:28pm UTC](https://discourse.julialang.org/t/slow-model-construction-in-jump/95455/15 "2023-03-14T22:28:03Z")

</div>

> [@odow](#):
>
> 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.

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`.