From e6d39f0650891b2d56e70be0944f02b2b1b41a5b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Mon, 21 Oct 2019 14:48:55 +0200 Subject: [PATCH 1/4] Simplify thanks for new API in CSDP --- deps/constants.jl | 4 +- src/CSDP.jl | 2 +- src/MOI_wrapper.jl | 102 +++++++++++++++++++++++++++++------------- src/blockmat.h.jl | 4 +- src/blockmat.jl | 15 +++---- src/declarations.h.jl | 63 +++++++++++++++++++++++--- src/declarations.jl | 79 +++++++++++++++++++++++++++++++- test/MOI_wrapper.jl | 6 ++- 8 files changed, 222 insertions(+), 53 deletions(-) diff --git a/deps/constants.jl b/deps/constants.jl index edaf66c..a55ed18 100644 --- a/deps/constants.jl +++ b/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") diff --git a/src/CSDP.jl b/src/CSDP.jl index c11ebe4..6aede6e 100644 --- a/src/CSDP.jl +++ b/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") diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 1b643be..747032d 100644 --- a/src/MOI_wrapper.jl +++ b/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,12 +82,14 @@ 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) @@ -91,12 +97,22 @@ function MOI.empty!(optimizer::Optimizer) 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,52 +211,78 @@ 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{blockmatrix}(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 @@ -248,9 +290,9 @@ function MOI.optimize!(optimizer::Optimizer) 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 diff --git a/src/blockmat.h.jl b/src/blockmat.h.jl index 415640b..e814c1c 100644 --- a/src/blockmat.h.jl +++ b/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 diff --git a/src/blockmat.jl b/src/blockmat.jl index e9b4ce2..3ca4378 100644 --- a/src/blockmat.jl +++ b/src/blockmat.jl @@ -15,7 +15,6 @@ end export fptr, ptr function mywrap(X::blockmatrix) - finalizer(free_blockmatrix, X) BlockMatrix(X) end @@ -28,7 +27,7 @@ function mywrap(x::Ptr{T}, len) where T # because the pointer it has has an offset y = _unsafe_wrap(Array, x + sizeof(T), len, false) # fptr takes care of this offset - finalizer(s -> Libc.free(fptr(s)), y) + #finalizer(s -> Libc.free(fptr(s)), y) y end @@ -79,10 +78,11 @@ function blockreczeros(n) end end -function Base.size(A::BlockRec) - n = A.csdp.blocksize +function Base.size(A::blockrec) + n = A.blocksize (n, n) end +Base.size(A::BlockRec) = size(A.csdp) function Base.getindex(A::BlockRec, i, j) n = A.csdp.blocksize if A.csdp.blockcategory == MATRIX @@ -309,16 +309,13 @@ function Base.getindex(A::Union{BlockMatrix, ConstraintMatrix}, i::Integer) block(A, i) end +block(A::blockmatrix, i::Integer) = getblockrec(A, i) function block(A::Union{BlockMatrix, ConstraintMatrix}, i::Integer) A.jblocks[i] end +nblocks(A::blockmatrix) = A.nblocks nblocks(A::Union{BlockMatrix, ConstraintMatrix}) = length(A.jblocks) -function free_blockmatrix(m::blockmatrix) - ccall((:free_mat, CSDP.csdp), Nothing, (blockmatrix,), m) -end -export free_blockmatrix - export BlockMatrix, ConstraintMatrix """Solver status""" diff --git a/src/declarations.h.jl b/src/declarations.h.jl index e78a015..5d914ec 100644 --- a/src/declarations.h.jl +++ b/src/declarations.h.jl @@ -179,6 +179,11 @@ function initparams(params::Ptr{paramstruc},pprintlevel::Ptr{Cint}) ccall((:initparams,CSDP.csdp),Nothing,(Ptr{paramstruc},Ptr{Cint}),params,pprintlevel) end +function loaded_initsoln(problem::Ptr{Cvoid}, X::Ref{blockmatrix}, Z::Ref{blockmatrix}) + y = Ref{Ptr{Cdouble}}(C_NULL) + ccall((:loaded_initsoln,CSDP.csdp),Nothing,(Ptr{Cvoid},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix}),problem,X,y,Z) + return y[] +end function initsoln(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix}) X = Ref{blockmatrix}(blockmatrix(0, C_NULL)) y = Ref{Ptr{Cdouble}}(C_NULL) @@ -260,6 +265,38 @@ function free_prob(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::P ccall((:free_prob,CSDP.csdp),Nothing,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},blockmatrix),n,k,C,a,constraints,X,y,Z) end +#function new_blockmatrix(nblocks::Integer) +# return ccall((:new_blockmatrix,CSDP.csdp),Ptr{Cvoid},(Cint,),nblocks) +#end +#function free_blockmatrix(C::Ptr{Cvoid}) +# ccall((:free_blockmatrix,CSDP.csdp),Nothing,(Ptr{Cvoid},),C) +#end + + +function Base.getindex(A::blockrec, i::Integer, j::Integer) + return ccall((:getindex,CSDP.csdp), Cdouble, (blockrec,Cint,Cint), A, i, j) +end +function getblockrec(A::blockmatrix, i::Integer) + return ccall((:getblockrec,CSDP.csdp), blockrec, (blockmatrix,Cint), A, i) +end + +function allocate_loading_prob(pC::Ref{blockmatrix}, block_dims::Ptr{Cint}, num_constraints::Integer, num_entries::Ptr{Cint}, printlevel::Integer) + return ccall((:allocate_loading_prob,CSDP.csdp),Ptr{Cvoid},(Ref{blockmatrix},Ptr{Cint},Cint,Ptr{Cint},Cint),pC,block_dims,num_constraints,num_entries,printlevel) +end +function free_loading_prob(problem::Ptr{Cvoid}) + ccall((:free_loading_prob,CSDP.csdp),Nothing,(Ptr{Cvoid},),problem) +end +function free_loaded_prob(problem::Ptr{Cvoid},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix) + ccall((:free_loaded_prob,CSDP.csdp),Nothing,(Ptr{Cvoid},blockmatrix,Ptr{Cdouble},blockmatrix),problem,X,y,Z) +end + +function setconstant(problem::Ptr{Cvoid}, mat::Integer, ent::Cdouble) + ccall((:setconstant,CSDP.csdp),Nothing,(Ptr{Cvoid},Cint,Cdouble),problem,mat,ent) +end +function addentry(problem::Ptr{Cvoid}, mat::Integer, blk::Integer, indexi::Integer, indexj::Integer, ent::Cdouble, allow_duplicates::Integer) + ccall((:addentry,CSDP.csdp),Cint,(Ptr{Cvoid},Cint,Cint,Cint,Cint,Cdouble,Cint),problem,mat,blk,indexi,indexj,ent,allow_duplicates) +end + function sdp(n::Cint, k::Cint, C::blockmatrix, a::Ptr{Cdouble}, constant_offset::Cdouble, constraints::Ptr{constraintmatrix}, byblocks::Ptr{Ptr{sparseblock}}, @@ -270,8 +307,8 @@ function sdp(n::Cint, k::Cint, O::Ptr{Cdouble}, rhs::Ptr{Cdouble}, dy::Ptr{Cdouble}, dy1::Ptr{Cdouble}, Fp::Ptr{Cdouble}, printlevel::Cint, parameters::paramstruc) - pobj = Ref{Cdouble}(.0) - dobj = Ref{Cdouble}(.0) + pobj = Ref{Cdouble}(0.0) + dobj = Ref{Cdouble}(0.0) work1 = blockmatrix(); alloc_mat(C, Ref{blockmatrix}(work1)); work2 = blockmatrix(); alloc_mat(C, Ref{blockmatrix}(work2)); @@ -333,14 +370,28 @@ function sdp(n::Cint, k::Cint, free_mat(dZ); free_mat(dX); - status, pobj[], dobj[] + return status, pobj[], dobj[] end +function loaded_sdp(problem::Ptr{Cvoid},constant_offset::Cdouble,pX::Ref{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ref{blockmatrix},printlevel::Cint,parameters::paramstruc) + pobj = Ref{Cdouble}(0.0) + dobj = Ref{Cdouble}(0.0) + status = ccall((:loaded_sdp,CSDP.csdp),Cint, + (Ptr{Cvoid},Cdouble, Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix},Ref{Cdouble},Ref{Cdouble},Cint,paramstruc), + problem, constant_offset,pX, py, pZ,pobj,dobj,printlevel,parameters) + return status, pobj[], dobj[] +end +function parametrized_sdp(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},constant_offset::Cdouble,pX::Ptr{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ptr{blockmatrix},printlevel::Cint,parameters::paramstruc) + pobj = Ref{Cdouble}(0.0) + dobj = Ref{Cdouble}(0.0) + status = ccall((:parametrized_sdp,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,Ptr{blockmatrix},Ref{Ptr{Cdouble}},Ptr{blockmatrix},Ref{Cdouble},Ref{Cdouble},Cint,paramstruc),n,k,C,a,constraints,constant_offset,pX,py,pZ,pobj,dobj,printlevel,parameters) + return status, pobj[], dobj[] +end function easy_sdp(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},constant_offset::Cdouble,pX::Ptr{blockmatrix},py::Ref{Ptr{Cdouble}},pZ::Ptr{blockmatrix}) - pobj = Ref{Cdouble}(.0) - dobj = Ref{Cdouble}(.0) + pobj = Ref{Cdouble}(0.0) + dobj = Ref{Cdouble}(0.0) status = ccall((:easy_sdp,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,Ptr{blockmatrix},Ref{Ptr{Cdouble}},Ptr{blockmatrix},Ref{Cdouble},Ref{Cdouble}),n,k,C,a,constraints,constant_offset,pX,py,pZ,pobj,dobj) - status, pobj[], dobj[] + return status, pobj[], dobj[] end function tweakgap(n::Cint,k::Cint,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},gap::Cdouble,Z::blockmatrix,dZ::blockmatrix,y::Ptr{Cdouble},dy::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,work4::blockmatrix,workvec1::Ptr{Cdouble},workvec2::Ptr{Cdouble},workvec3::Ptr{Cdouble},workvec4::Ptr{Cdouble},printlevel::Cint) diff --git a/src/declarations.jl b/src/declarations.jl index f335ede..06308df 100644 --- a/src/declarations.jl +++ b/src/declarations.jl @@ -1,5 +1,31 @@ export initsoln, easy_sdp +struct LoadingProblem + ptr::Ptr{Cvoid} +end + +function allocate_loading_prob(pC::Ref{blockmatrix}, block_dims::Vector{Cint}, num_constraints::Integer, num_entries::Matrix{Cint}, printlevel::Integer) + ptr = allocate_loading_prob(pC, fptr(block_dims), num_constraints, pointer(num_entries), printlevel) + return LoadingProblem(ptr) +end +free_loading_prob(problem::LoadingProblem) = free_loading_prob(problem.ptr) +function free_loaded_prob(problem::LoadingProblem, X::blockmatrix, y::Vector{Cdouble}, Z::blockmatrix) + free_loaded_prob(problem.ptr, X, fptr(y), Z) +end + +function setconstant(problem::LoadingProblem, mat::Integer, ent::Cdouble) + setconstant(problem.ptr, mat, ent) +end + +function addentry(problem::LoadingProblem, mat::Integer, blk::Integer, indexi::Integer, indexj::Integer, ent::Cdouble, allow_duplicates::Bool) + ret = addentry(problem.ptr, mat, blk, indexi, indexj, ent, allow_duplicates) + return !iszero(ret) +end + +function loaded_initsoln(problem::LoadingProblem, num_constraints::Integer, X::Ref{blockmatrix}, Z::Ref{blockmatrix}) + y = loaded_initsoln(problem.ptr, X, Z) + return mywrap(y, num_constraints) +end function initsoln(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}) m = length(As) X, y, Z = initsoln(Cint(size(C, 1)), Cint(m), C.csdp, fptr(b), fptr(As)) @@ -90,9 +116,58 @@ function sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X @assert Xcsdp == X.csdp @assert ycsdp == fptr(y) @assert Zcsdp == Z.csdp - status, pobj, dobj + return status, pobj, dobj +end + +function loaded_sdp(problem::LoadingProblem, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, params::Dict{Symbol}) + return loaded_sdp(problem, X, y, Z, options(params)...) +end +function loaded_sdp(problem::LoadingProblem, X::Ref{blockmatrix}, y::Vector{Cdouble}, Z::Ref{blockmatrix}, printlevel, params) + # I pass pointers py to X, y and Z but only *pX, *py and *pZ are + # used in the code so no need to worry, they won't change :) + ycsdp = Ref{Ptr{Cdouble}}(fptr(y)) + + status, pobj, dobj = loaded_sdp( + problem.ptr, # problem::Ptr{Cvoid} + 0.0, # constant_offset::Cdouble + X, # pX::Ptr{blockmatrix} + ycsdp, # py::Ptr{Cdouble} + Z, # pZ::Ptr{blockmatrix} + Cint(printlevel), # printlevel::Cint + params) # parameters::paramstruc + # I can even assert that they won't change + @assert ycsdp[] == fptr(y) + return status, pobj, dobj +end +function parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, params::Dict{Symbol}) + return parametrized_sdp(C, b, As, X, y, Z, options(params)...) end +function parametrized_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{ConstraintMatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix, printlevel, params) + # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are + # used in the code so no need to worry, they won't change :) + Xcsdp = X.csdp + ycsdp = Ref{Ptr{Cdouble}}(fptr(y)) + Zcsdp = Z.csdp + Ascsdp = map(A->A.csdp, As) + status, pobj, dobj = parametrized_sdp( + Cint(size(C, 1)), # n::Cint + Cint(length(As)), # k::Cint + C.csdp, # C::blockmatrix + fptr(b), # a::Ptr{Cdouble} + fptr(Ascsdp), # constraints::Ptr{constraintmatrix} + 0.0, # constant_offset::Cdouble + ptr(Xcsdp), # pX::Ptr{blockmatrix} + ycsdp, # py::Ptr{Cdouble} + ptr(Zcsdp), # pZ::Ptr{blockmatrix} + Cint(printlevel), # printlevel::Cint + params) # parameters::paramstruc + # I can even assert that they won't change + @assert Xcsdp == X.csdp + @assert ycsdp[] == fptr(y) + @assert Zcsdp == Z.csdp + return status, pobj, dobj +end function easy_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}, X::BlockMatrix, y::Vector{Cdouble}, Z::BlockMatrix) # I pass pointers pX, py and pZ to X, y and Z but only *pX, *py and *pZ are @@ -114,7 +189,7 @@ function easy_sdp(C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatri @assert Xcsdp == X.csdp @assert ycsdp[] == fptr(y) @assert Zcsdp == Z.csdp - status, pobj, dobj + return status, pobj, dobj end function write_prob(fname::String, C::BlockMatrix, b::Vector{Cdouble}, As::Vector{constraintmatrix}) diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 973e45c..6e61d4e 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -54,13 +54,17 @@ end MOIB.remove_bridge(bridged, MOIB.Constraint.ScalarSlackBridge{Float64}) MOIT.contlineartest(bridged, config, [ # Finds `MOI.ALMOST_OPTIMAL` instead of `MOI.OPTIMAL` - "linear10b" + "linear10b", + # Empty constraint + "linear15" ]) end @testset "Continuous Conic" begin MOIT.contconictest(bridged, config, [ # Finds `MOI.OPTIMAL` instead of `MOI.INFEASIBLE`. "soc3", + # Empty constraint `c4` + "psdt2", # See https://github.com/coin-or/Csdp/issues/11 "rotatedsoc1v", # Missing bridges From 4fce9324640b6ff9519478b0d6685ecad35a164b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 22 Oct 2019 16:01:18 +0200 Subject: [PATCH 2/4] Update MPB wrapper to the new CSDP interface --- src/MOI_wrapper.jl | 2 +- src/MPB_wrapper.jl | 74 ++++++++++++++++++++++++++++++------------- src/blockmat.jl | 2 +- src/declarations.h.jl | 10 +++++- 4 files changed, 63 insertions(+), 25 deletions(-) diff --git a/src/MOI_wrapper.jl b/src/MOI_wrapper.jl index 747032d..c9ebd68 100644 --- a/src/MOI_wrapper.jl +++ b/src/MOI_wrapper.jl @@ -219,7 +219,7 @@ function MOIU.load_variables(optimizer::Optimizer, nvars) for (key, value) in optimizer.num_entries num_entries[key...] = value end - optimizer.problem = allocate_loading_prob(Ref{blockmatrix}(optimizer.C), optimizer.blockdims, length(optimizer.b), num_entries, 3) + 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 duplicate = addentry(optimizer.problem, 1, length(optimizer.blockdims), 1, 1, 1.0, true) diff --git a/src/MPB_wrapper.jl b/src/MPB_wrapper.jl index 64caa6a..a268dde 100644 --- a/src/MPB_wrapper.jl +++ b/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) diff --git a/src/blockmat.jl b/src/blockmat.jl index 3ca4378..a31545c 100644 --- a/src/blockmat.jl +++ b/src/blockmat.jl @@ -305,7 +305,7 @@ function ConstraintMatrix(csdp::constraintmatrix, k::Integer) end # Needed by MPB_wrapper -function Base.getindex(A::Union{BlockMatrix, ConstraintMatrix}, i::Integer) +function Base.getindex(A::Union{blockmatrix, BlockMatrix, ConstraintMatrix}, i::Integer) block(A, i) end diff --git a/src/declarations.h.jl b/src/declarations.h.jl index 5d914ec..26290ae 100644 --- a/src/declarations.h.jl +++ b/src/declarations.h.jl @@ -216,7 +216,15 @@ function read_sol(fname::Ptr{UInt8},n::Cint,k::Cint,C::blockmatrix,pX::Ptr{block ccall((:read_sol,CSDP.csdp),Cint,(Ptr{UInt8},Cint,Cint,blockmatrix,Ptr{blockmatrix},Ptr{Ptr{Cdouble}},Ptr{blockmatrix}),fname,n,k,C,pX,py,pZ) end -function read_prob(fname::String,printlevel::Integer=0) +function load_prob_from_file(fname::String,C::Ref{blockmatrix},printlevel::Integer=1) + problem = Ref{Ptr{Cvoid}}(C_NULL) + ret = ccall((:load_prob_from_file,CSDP.csdp),Cint,(Ptr{UInt8},Ref{blockmatrix},Ref{Ptr{Cvoid}},Cint),fname,C,problem,printlevel) + if !iszero(ret) + error("`CSDP.load_prob_from_file` failed.") + end + return LoadingProblem(problem[]) +end +function read_prob(fname::String,printlevel::Integer=1) n = Ref{Cint}(0) k = Ref{Cint}(0) C = Ref{blockmatrix}(blockmatrix(0, C_NULL)) From 97928e83943b9a6a61c05013d69bbf0a8736af73 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Tue, 22 Oct 2019 22:01:24 +0200 Subject: [PATCH 3/4] Exclude solve_unbounded_model --- test/MOI_wrapper.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/MOI_wrapper.jl b/test/MOI_wrapper.jl index 6e61d4e..4c8b981 100644 --- a/test/MOI_wrapper.jl +++ b/test/MOI_wrapper.jl @@ -33,6 +33,8 @@ end @testset "Unit" begin MOIT.unittest(bridged, config, [ + # `NUMERICAL_ERROR` on Mac: https://travis-ci.org/JuliaOpt/CSDP.jl/jobs/601302777#L217-L219 + "solve_unbounded_model", # `NumberOfThreads` not supported. "number_threads", # `TimeLimitSec` not supported. From 8f63bbf44f52ffff4463dd65ab211454ce693f3b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Beno=C3=AEt=20Legat?= Date: Wed, 23 Oct 2019 11:05:10 +0200 Subject: [PATCH 4/4] Exclude MPB SOC tests --- test/MPB_wrapper.jl | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/test/MPB_wrapper.jl b/test/MPB_wrapper.jl index d1249f4..4948c8a 100644 --- a/test/MPB_wrapper.jl +++ b/test/MPB_wrapper.jl @@ -17,9 +17,10 @@ end coniclineartest(solver, duals=true, tol=1e-6) end -@testset "Conic SOC tests" begin - conicSOCtest(CSDP.CSDPSolver(printlevel=0, write_prob="soc.prob"), duals=true, tol=1e-6) -end +# :Error for SOC1 on Mac OS: https://travis-ci.org/JuliaOpt/CSDP.jl/jobs/601461403#L389-L391 +#@testset "Conic SOC tests" begin +# conicSOCtest(CSDP.CSDPSolver(printlevel=0, write_prob="soc.prob"), duals=true, tol=1e-6) +#end # CSDP returns :Suboptimal for SOCRotated1 # @testset "Conic SOC rotated tests" begin