Merge 4fce932 into e367784

deps/constants.jl

@@ -3,9 +3,9 @@ const JULIA_LAPACK = false
 const suffix       = JULIA_LAPACK ? ".64" : ""
 const version      = "6.2.0"
 const libname      = "libcsdp$suffix.$(Libdl.dlext)"
-const csdpversion  = "Csdp-$version"
+const csdpversion  = "Csdp-readprob"
 const download_url =
-    "http://www.coin-or.org/download/source/Csdp/Csdp-$version.tgz"
+    "https://github.com/blegat/Csdp/archive/readprob.zip"
 
 patchf      = joinpath(dirname(@__FILE__), "src$suffix", "debug-mat.c")
 srcdir      = joinpath(dirname(@__FILE__), "src$suffix", csdpversion, "lib")

src/CSDP.jl

@@ -12,8 +12,8 @@ end
 
 export Blockmatrix
 
-include("blockmat.h.jl")
 include("blockdiag.jl")
+include("blockmat.h.jl")
 include("blockmat.jl")
 include("declarations.h.jl")
 include("declarations.jl")

src/MOI_wrapper.jl

@@ -6,14 +6,15 @@ const AFFEQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo
 mutable struct Optimizer <: MOI.AbstractOptimizer
     objconstant::Cdouble
     objsign::Int
-    blockdims::Vector{Int}
+    blockdims::Vector{Cint}
     varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
+    num_entries::Dict{Tuple{Int, Int}, Int}
     b::Vector{Cdouble}
-    C::Union{Nothing, BlockMatrix}
-    As::Union{Nothing, Vector{ConstraintMatrix}}
-    X::Union{Nothing, BlockMatrix}
+    C::blockmatrix
+    problem::Union{Nothing, LoadingProblem}
+    X::blockmatrix
     y::Union{Nothing, Vector{Cdouble}}
-    Z::Union{Nothing, BlockMatrix}
+    Z::blockmatrix
     status::Cint
     pobj::Cdouble
     dobj::Cdouble
@@ -22,12 +23,15 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     options::Dict{Symbol, Any}
     function Optimizer(; kwargs...)
         optimizer = new(
-            zero(Cdouble), 1, Int[], Tuple{Int, Int, Int}[], Cdouble[],
-            nothing, nothing, nothing, nothing, nothing,
+            zero(Cdouble), 1, Cint[], Tuple{Int, Int, Int}[],
+            Dict{Tuple{Int, Int}, Int}(), Cdouble[],
+            blockmatrix(), nothing, blockmatrix(), nothing, blockmatrix(),
             -1, NaN, NaN, NaN, false, Dict{Symbol, Any}())
         for (key, value) in kwargs
             MOI.set(optimizer, MOI.RawParameter(key), value)
         end
+        # May need to call `free_loaded_prob` and `free_loading_prob`.
+        finalizer(MOI.empty!, optimizer)
         return optimizer
     end
 end
@@ -78,25 +82,37 @@ end
 
 function MOI.is_empty(optimizer::Optimizer)
     return iszero(optimizer.objconstant) &&
-        optimizer.objsign == 1 &&
+        isone(optimizer.objsign) &&
         isempty(optimizer.blockdims) &&
         isempty(optimizer.varmap) &&
+        isempty(optimizer.num_entries) &&
         isempty(optimizer.b) &&
-        optimizer.C === nothing &&
-        optimizer.As === nothing
+        iszero(optimizer.C.nblocks) &&
+        optimizer.C.blocks == C_NULL &&
+        optimizer.problem === nothing
 end
 
 function MOI.empty!(optimizer::Optimizer)
     optimizer.objconstant = zero(Cdouble)
     optimizer.objsign = 1
     empty!(optimizer.blockdims)
     empty!(optimizer.varmap)
+    empty!(optimizer.num_entries)
     empty!(optimizer.b)
-    optimizer.C = nothing
-    optimizer.As = nothing
-    optimizer.X = nothing
+    if optimizer.problem !== nothing
+        if optimizer.y !== nothing
+            free_loaded_prob(optimizer.problem, optimizer.X, optimizer.y, optimizer.Z)
+        end
+        free_loading_prob(optimizer.problem)
+    end
+    optimizer.problem = nothing
+    optimizer.C.nblocks = 0
+    optimizer.C.blocks = C_NULL
+    optimizer.X.nblocks = 0
+    optimizer.X.blocks = C_NULL
     optimizer.y = nothing
-    optimizer.Z = nothing
+    optimizer.Z.nblocks = 0
+    optimizer.Z.blocks = C_NULL
     optimizer.status = -1
     optimizer.pobj = 0.0
     optimizer.dobj = 0.0
@@ -140,12 +156,12 @@ function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) e
 # Loads objective coefficient α * vi
 function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex)
     blk, i, j = varmap(optimizer, vi)
+    # in SDP format, it is max and in MPB Conic format it is min
     coef = optimizer.objsign * α
     if i != j
         coef /= 2
     end
-    # in SDP format, it is max and in MPB Conic format it is min
-    block(optimizer.C, blk)[i, j] = coef
+    addentry(optimizer.problem, 0, blk, i, j, coef, true)
 end
 function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
     obj = MOIU.canonical(f)
@@ -195,62 +211,88 @@ function MOIU.load_variables(optimizer::Optimizer, nvars)
     if dummy
         # See https://github.com/coin-or/Csdp/issues/2
         optimizer.b = [one(Cdouble)]
-        optimizer.blockdims = [optimizer.blockdims; -1]
+        optimizer.blockdims = [optimizer.blockdims; Cint(-1)]
+        count_entry(optimizer, 1, length(optimizer.blockdims))
     end
-    optimizer.C = blockmatzeros(optimizer.blockdims)
-    optimizer.As = [constrmatzeros(i, optimizer.blockdims) for i in eachindex(optimizer.b)]
+    optimizer.C.nblocks = length(optimizer.blockdims)
+    num_entries = zeros(Cint, length(optimizer.b), length(optimizer.blockdims))
+    for (key, value) in optimizer.num_entries
+        num_entries[key...] = value
+    end
+    optimizer.problem = allocate_loading_prob(Ref(optimizer.C), optimizer.blockdims, length(optimizer.b), num_entries, 3)
     if dummy
         # See https://github.com/coin-or/Csdp/issues/2
-        block(optimizer.As[1], length(optimizer.blockdims))[1, 1] = 1
+        duplicate = addentry(optimizer.problem, 1, length(optimizer.blockdims), 1, 1, 1.0, true)
+        @assert !duplicate
     end
 
 end
 
+function count_entry(optimizer::Optimizer, con_idx::Integer, blk::Integer)
+    key = (con_idx, blk)
+    optimizer.num_entries[key] = get(optimizer.num_entries, key, 0) + 1
+end
+
 function MOIU.allocate_constraint(optimizer::Optimizer,
                                   func::MOI.ScalarAffineFunction{Cdouble},
                                   set::MOI.EqualTo{Cdouble})
+    if !iszero(MOI.constant(func))
+        throw(MOI.ScalarFunctionConstantNotZero{
+            Cdouble, MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}(
+                MOI.constant(func)))
+    end
     push!(optimizer.b, MOI.constant(set))
+    func = MOIU.canonical(func) # sum terms with same variables and same output_index
+    for t in func.terms
+        if !iszero(t.coefficient)
+            blk, i, j = varmap(optimizer, t.variable_index)
+            count_entry(optimizer, length(optimizer.b), blk)
+        end
+    end
     return AFFEQ(length(optimizer.b))
 end
 
-function MOIU.load_constraint(m::Optimizer, ci::AFFEQ,
+function MOIU.load_constraint(optimizer::Optimizer, ci::AFFEQ,
                               f::MOI.ScalarAffineFunction, s::MOI.EqualTo)
     if !iszero(MOI.constant(f))
         throw(MOI.ScalarFunctionConstantNotZero{
             Cdouble, MOI.ScalarAffineFunction{Cdouble}, MOI.EqualTo{Cdouble}}(
                 MOI.constant(f)))
     end
-    f = MOIU.canonical(f) # sum terms with same variables and same outputindex
+    setconstant(optimizer.problem, ci.value, MOI.constant(s))
+    f = MOIU.canonical(f) # sum terms with same variables and same output_index
+    if isempty(f.terms)
+        throw(ArgumentError("Empty constraint $ci: $f-in-$s. Not supported by CSDP."))
+    end
     for t in f.terms
         if !iszero(t.coefficient)
-            blk, i, j = varmap(m, t.variable_index)
+            blk, i, j = varmap(optimizer, t.variable_index)
             coef = t.coefficient
             if i != j
                 coef /= 2
             end
-            block(m.As[ci.value], blk)[i, j] = coef
+            duplicate = addentry(optimizer.problem, ci.value, blk, i, j, coef, true)
+            @assert !duplicate
         end
     end
 end
 
 
 function MOI.optimize!(optimizer::Optimizer)
-    As = map(A->A.csdp, optimizer.As)
-
     write_prob(optimizer)
 
     start_time = time()
-    optimizer.X, optimizer.y, optimizer.Z = initsoln(optimizer.C, optimizer.b, As)
+    optimizer.y = loaded_initsoln(optimizer.problem, length(optimizer.b), Ref(optimizer.X), Ref(optimizer.Z))
 
     options = optimizer.options
     if optimizer.silent
         options = copy(options)
         options[:printlevel] = 0
     end
 
-    optimizer.status, optimizer.pobj, optimizer.dobj = sdp(
-        optimizer.C, optimizer.b, optimizer.As, optimizer.X, optimizer.y,
-        optimizer.Z, options)
+    optimizer.status, optimizer.pobj, optimizer.dobj = loaded_sdp(
+        optimizer.problem, Ref(optimizer.X), optimizer.y,
+        Ref(optimizer.Z), options)
     optimizer.solve_time = time() - start_time
 end
 

src/MPB_wrapper.jl

@@ -20,18 +20,21 @@ end
 CSDPSolver(; kwargs...) = CSDPSolver(checkoptions(Dict{Symbol,Any}(kwargs)))
 
 mutable struct CSDPMathProgModel <: SDM.AbstractSDModel
-    C
-    b
-    As
-    X
-    y
-    Z
+    blockdims::Vector{Cint}
+    b::Vector{Cdouble}
+    entries::Vector{Tuple{Cint,Cint,Cint,Cint,Cdouble}}
+    C::blockmatrix
+    problem::Union{Nothing, LoadingProblem}
+    X::blockmatrix
+    y::Union{Nothing, Vector{Cdouble}}
+    Z::blockmatrix
     status::Cint
     pobj::Cdouble
     dobj::Cdouble
     options::Dict{Symbol,Any}
     function CSDPMathProgModel(; kwargs...)
-        new(nothing, nothing, nothing, nothing, nothing, nothing,
+        new(Cint[], Cdouble[], Tuple{Cint,Cint,Cint,Cint,Cdouble}[],
+            blockmatrix(), nothing, blockmatrix(), nothing, blockmatrix(),
             -1, 0.0, 0.0, checkoptions(Dict{Symbol, Any}(kwargs)))
     end
 end
@@ -47,39 +50,66 @@ function MPB.setvartype!(m::CSDPMathProgModel, vtype, blk, i, j)
 end
 
 function MPB.loadproblem!(m::CSDPMathProgModel, filename::AbstractString)
+    if m.problem !== nothing
+        if m.y !== nothing
+            free_loaded_prob(m.problem, m.X, m.y, m.Z)
+        end
+        free_loading_prob(m.problem)
+    end
+    m.problem = nothing
+    m.y = nothing
     if endswith(filename,".dat-s")
-       m.C, m.b, As = read_prob(filename)
-       m.As = [ConstraintMatrix(As[i], i) for i in 1:length(As)]
+        m.problem = load_prob_from_file(filename, Ref(m.C))
     else
-       error("unrecognized input format extension in $filename")
+        error("unrecognized input format extension in $filename")
     end
 end
 #writeproblem(m, filename::String)
 function MPB.loadproblem!(m::CSDPMathProgModel, blkdims::Vector{Int}, constr::Int)
-    m.C = blockmatzeros(blkdims)
+    if m.problem !== nothing
+        if m.y !== nothing
+            free_loaded_prob(m.problem, m.X, m.y, m.Z)
+        end
+        free_loading_prob(m.problem)
+    end
+    m.problem = nothing
+    m.y = nothing
+    m.blockdims = blkdims
     m.b = zeros(Cdouble, constr)
-    m.As = [constrmatzeros(i, blkdims) for i in 1:constr]
+    empty!(m.entries)
 end
 
 function SDM.setconstrB!(m::CSDPMathProgModel, val, constr::Integer)
     m.b[constr] = val
 end
 function SDM.setconstrentry!(m::CSDPMathProgModel, coef, constr::Integer, blk::Integer, i::Integer, j::Integer)
-    block(m.As[constr], blk)[i, j] = coef
+    push!(m.entries, (constr, blk, i, j, coef))
 end
 function SDM.setobjentry!(m::CSDPMathProgModel, coef, blk::Integer, i::Integer, j::Integer)
-    block(m.C, blk)[i, j] = coef
+    push!(m.entries, (0, blk, i, j, coef))
 end
 
 function MPB.optimize!(m::CSDPMathProgModel)
-    As = map(A->A.csdp, m.As)
-
-    write_prob(m)
-
-    m.X, m.y, m.Z = initsoln(m.C, m.b, As)
-    #verbose = get(m.options, :verbose, true)
-    #m.status, m.pobj, m.dobj = easy_sdp(m.C, m.b, As, m.X, m.y, m.Z, verbose)
-    m.status, m.pobj, m.dobj = sdp(m.C, m.b, m.As, m.X, m.y, m.Z, m.options)
+    if m.problem === nothing
+        # `m.problem` is not `nothing` if it was loaded from a file.
+        m.C.nblocks = length(m.blockdims)
+        num_entries = zeros(Cint, length(m.b), length(m.blockdims))
+        for entry in m.entries
+            if entry[1] > 0
+                num_entries[entry[1], entry[2]] += 1
+            end
+        end
+        m.problem = allocate_loading_prob(Ref(m.C), m.blockdims, length(m.b), num_entries, 3)
+        for (i, x) in enumerate(m.b)
+            setconstant(m.problem, i, x)
+        end
+        for entry in m.entries
+            duplicate = addentry(m.problem, entry..., true)
+            @assert !duplicate
+        end
+    end
+    m.y = loaded_initsoln(m.problem, length(m.b), Ref(m.X), Ref(m.Z))
+    m.status, m.pobj, m.dobj = loaded_sdp(m.problem, Ref(m.X), m.y, Ref(m.Z), m.options)
 end
 
 function MPB.status(m::CSDPMathProgModel)

src/blockmat.h.jl

@@ -20,13 +20,13 @@ struct blockdatarec
     _blockdatarec::Ptr{Cdouble}
 end
 
-struct blockrec
+struct blockrec <: AbstractMatrix{Cdouble}
     data::blockdatarec
     blockcategory::blockcat
     blocksize::csdpshort
 end
 
-mutable struct blockmatrix
+mutable struct blockmatrix <: AbstractBlockMatrix{Cdouble}
     nblocks::Cint
     blocks::Ptr{blockrec}
 end