Construct LP from matrix using GRBnewmodel, GRBaddvars, and GRBaddconstr

#=

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

returnvalue.(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))