Demo code for Stéphan Laurent blog post comment

glpkblog1.jl

import GLPK # NOTE: it's better to use 'import' with GLPK,
            #       rather than 'using'

μ = [1/7, 2/7, 4/7]
ν = [1/4, 1/4, 1/2]

n = length(μ)
D = 1 .- eye(n) # NOTE: starting from julia v0.3, you won't be
                #       allowed to use '-' here, you need to
                #       use '.-'

# alternative definition, only works in julia v0.3:
# D = ones(n,n) - I

M1 = zeros(n, n^2)
for i = 1:n
    M1[i, (1:n)+n*(i-1)] = 1
end
M2 = repmat(eye(n), 1, n)
A = vcat(M1, M2)

lp = GLPK.Prob()
GLPK.set_prob_name(lp, "kanto")

GLPK.set_obj_dir(lp, GLPK.MIN)

GLPK.add_rows(lp, 2n)
for i = 1:n
    GLPK.set_row_bnds(lp, i, GLPK.FX, μ[i], μ[i])
    GLPK.set_row_bnds(lp, i+n, GLPK.FX, ν[i], ν[i])
end

GLPK.add_cols(lp, n^2)
k = 0
for i = 1:n, j = 1:n # NOTE: nested loops can be written like this
    k += 1
    GLPK.set_col_bnds(lp, k, GLPK.LO, 0.0, 0.0)
    GLPK.set_obj_coef(lp, k, D[i,j])
end

GLPK.load_matrix(lp, sparse(A)) # NOTE: GLPK.load_matrix can work with sparse matrices
flag = GLPK.simplex(lp)
# flag = GLPK.exact(lp)
flag == 0 || error("GLPK failed")

result = GLPK.get_obj_val(lp)

glpkblog2.jl

using MathProgBase
# using GLPKMathProgInterface # NOTE: only needed to specify solver options
                              #       e.g. method = :Exact

μ = [1/7, 2/7, 4/7]
ν = [1/4, 1/4, 1/2]

n = length(μ)
D = 1 .- eye(n) # NOTE: starting from julia v0.3, you won't be
                #       allowed to use '-' here, you need to
                #       use '.-'

# alternative definition, only works in julia v0.3:
# D = ones(n,n) - I

M1 = zeros(n, n^2)
for i = 1:n
    M1[i, (1:n)+n*(i-1)] = 1
end
M2 = repmat(eye(n), 1, n)

A = vcat(M1, M2)
μν = vcat(μ, ν)
vD = vec(D)

sol = linprog(vD, A, '=', μν, 0.0, Inf)
#sol = linprog(vD, A, '=', μν, 0.0, Inf, GLPKSolverLP(method = :Exact))

sol.status == :Optimal || error("Optimization failed: $(sol.status)")

result = sol.objval