Merge pull request #63 from JuliaOpt/ml/memory

address memory corruption issue

Author: mlubin

src/ECOS.jl

@@ -43,9 +43,23 @@ function __init__()
     end
 end
 
-include("ECOSSolverInterface.jl")  # MathProgBase interface
 include("types.jl")  # All the types and constants defined in ecos.h
 
+# A julia container for matrices passed to ECOS.
+# This makes it easy for the Julia GC to track these arrays.
+# The format matches ECOS's spmat: 0-based compressed sparse column (CSC)
+struct ECOSMatrix
+    pr::Vector{Cdouble} # nzval
+    jc::Vector{Clong}   # colptr
+    ir::Vector{Clong}   # rowval
+end
+
+function ECOSMatrix(mat::AbstractMatrix)
+    sparsemat = sparse(mat)
+    return ECOSMatrix(sparsemat.nzval, sparsemat.colptr .- 1, sparsemat.rowval .- 1)
+end
+
+
 # setup  (direct interface)
 # Provide ECOS with a problem in the form
 # min c'x
@@ -66,35 +80,33 @@ include("types.jl")  # All the types and constants defined in ecos.h
 #               e.g. cone 1 is indices 4:6, cone 2 is indices 7:10
 #                    ->  q = [3, 4]
 #   e       Number of exponential cones present in problem
-#   Gpr, Gjc, Gir
-#           Non-zeros, column indices, and the row index arrays, respectively,
-#           for the matrix G represented in column compressed storage (CCS) format
-#           which is equivalent to Julia's SparseMatrixCSC
-#   Apr, Ajc, Air
-#           Equivalent to above for the matrix A.
+#   G       The G matrix in ECOSMatrix format.
+#   A       The A matrix in ECOSMatrix format.
 #           Can be all nothing if no equalities are present.
 #   c       Objective coefficients, length(c) == n
 #   h       RHS for inequality constraints, length(h) == m
 #   b       RHS for equality constraints, length(b) == b (can be nothing)
 # Returns a pointer to the ECOS pwork structure (Cpwork in ECOS.jl). See
 # types.jl for more information.
+# **NOTE**: ECOS retains references to the problem data passed in here.
+# You *must* ensure that G, A, c, h, and b are not freed until after cleanup(), otherwise
+# memory corruption will occur.
 function setup(n::Int, m::Int, p::Int, l::Int, ncones::Int, q::Union{Vector{Int},Void}, e::Int,
-                Gpr::Vector{Float64}, Gjc::Vector{Int}, Gir::Vector{Int},
-                Apr::Union{Vector{Float64},Void}, Ajc::Union{Vector{Int},Void}, Air::Union{Vector{Int},Void},
-                c::Vector{Float64}, h::Vector{Float64}, b::Union{Vector{Float64},Void}; kwargs...)
+               G::ECOSMatrix, A::Union{ECOSMatrix, Void},
+               c::Vector{Float64}, h::Vector{Float64}, b::Union{Vector{Float64},Void}; kwargs...)
     # Convert to canonical forms
     q = (q == nothing) ? convert(Ptr{Clong}, C_NULL) : convert(Vector{Clong},q)
-    Apr = (Apr == nothing) ? convert(Ptr{Cdouble}, C_NULL) : Apr
-    Ajc = (Ajc == nothing) ? convert(Ptr{Cdouble}, C_NULL) : convert(Vector{Clong},Ajc)
-    Air = (Air == nothing) ? convert(Ptr{Cdouble}, C_NULL) : convert(Vector{Clong},Air)
+    Apr = (A == nothing) ? convert(Ptr{Cdouble}, C_NULL) : A.pr
+    Ajc = (A == nothing) ? convert(Ptr{Cdouble}, C_NULL) : A.jc
+    Air = (A == nothing) ? convert(Ptr{Cdouble}, C_NULL) : A.ir
     b = (b == nothing) ? convert(Ptr{Cdouble}, C_NULL) : b
     problem_ptr = ccall((:ECOS_setup, ECOS.ecos), Ptr{Cpwork},
         (Clong, Clong, Clong, Clong, Clong, Ptr{Clong}, Clong,
          Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong},
          Ptr{Cdouble}, Ptr{Clong}, Ptr{Clong},
          Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}),
         n, m, p, l, ncones, q, e,
-        Gpr, convert(Vector{Clong},Gjc), convert(Vector{Clong},Gir),
+        G.pr, G.jc, G.ir,
         Apr, Ajc, Air,
         c, h, b)
 
@@ -122,47 +134,6 @@ function setup(n::Int, m::Int, p::Int, l::Int, ncones::Int, q::Union{Vector{Int}
     problem_ptr
 end
 
-# setup  (more general interface)
-# A more tolerant version of the above method that doesn't require
-# user to fiddle with internals of the sparse matrix format
-# User can pass nothing as argument for A, b, and q
-function setup(n, m, p, l, ncones, q, e, G, A, c, h, b; options...)
-    @assert m == l + sum(q) + 3e
-    @assert length(c) == n
-    @assert ncones == length(q)
-    if A == nothing
-        if b != nothing
-            @assert length(b) == 0
-            b = nothing
-        end
-        @assert p == 0
-        Apr = nothing
-        Ajc = nothing
-        Air = nothing
-    else
-        numrow, numcol = size(A)
-        @assert numcol == n
-        @assert numrow == length(b)
-        @assert numrow == p
-        sparseA = sparse(A)
-        Apr = convert(Vector{Float64}, sparseA.nzval)
-        Ajc = sparseA.colptr - 1  # C is 0-based
-        Air = sparseA.rowval - 1
-    end
-
-    @assert size(G) == (m, n)
-    @assert m == length(h)
-    sparseG = sparse(G)
-    Gpr = convert(Vector{Float64}, sparseG.nzval)
-    Gjc = sparseG.colptr - 1  # C is 0-based
-    Gir = sparseG.rowval - 1
-
-    setup(  n, m, p, l, ncones, q, e,
-            Gpr, Gjc, Gir,
-            Apr, Ajc, Air,
-            c, h, b; options...)
-end
-
 # solve
 # Solves the provided problem. Results are stored inside the structure,
 # but currently there is no convenient interface-provided way to access
@@ -178,4 +149,6 @@ function cleanup(problem::Ptr{Cpwork}, keepvars::Int = 0)
     ccall((:ECOS_cleanup, ECOS.ecos), Void, (Ptr{Cpwork}, Clong), problem, keepvars)
 end
 
+include("ECOSSolverInterface.jl")  # MathProgBase interface
+
 end # module

src/ECOSSolverInterface.jl

@@ -26,8 +26,8 @@ mutable struct ECOSMathProgModel <: AbstractConicModel
     ncones::Int                         # Number of SO cones
     conedims::Vector{Int}               # Dimension of each SO cone
     nexp_cones::Int                     # Number of exponential cones
-    G::SparseMatrixCSC{Float64,Int}     # The G matrix (inequalties)
-    A::SparseMatrixCSC{Float64,Int}     # The A matrix (equalities)
+    G::ECOSMatrix                       # The G matrix (inequalties)
+    A::ECOSMatrix                       # The A matrix (equalities)
     c::Vector{Float64}                  # The objective coeffs (always min)
     orig_sense::Symbol                  # Original objective sense
     h::Vector{Float64}                  # RHS for inequality
@@ -55,8 +55,8 @@ mutable struct ECOSMathProgModel <: AbstractConicModel
 end
 ECOSMathProgModel(;kwargs...) = ECOSMathProgModel(0,0,0,0,0,0,
                                         Int[],0,
-                                        spzeros(0,0),
-                                        spzeros(0,0),
+                                        ECOSMatrix(spzeros(0,0)),
+                                        ECOSMatrix(spzeros(0,0)),
                                         Float64[], :Min,
                                         Float64[], Float64[],
                                         :NotSolved, 0.0, 0.0,
@@ -85,8 +85,9 @@ function optimize!(m::ECOSMathProgModel)
         m.nvar, m.nineq, m.neq,
         m.npos, m.ncones, m.conedims, m.nexp_cones,
         m.G, m.A,
-        m.c[:],  # Seems to modify this
+        m.c,
         m.h, m.b; m.options...)
+    # Note: ECOS modifies problem data in setup() and restores it on cleanup()
 
     t = time()
     flag = solve(ecos_prob_ptr)
@@ -109,8 +110,9 @@ function optimize!(m::ECOSMathProgModel)
     m.primal_sol = unsafe_wrap(Array,ecos_prob.x, m.nvar)[:]
     m.dual_sol_eq   = unsafe_wrap(Array,ecos_prob.y, m.neq)[:]
     m.dual_sol_ineq = unsafe_wrap(Array,ecos_prob.z, m.nineq)[:]
-    m.obj_val = dot(m.c, m.primal_sol) * (m.orig_sense == :Max ? -1 : +1)
     cleanup(ecos_prob_ptr, 0)
+    # Compute this value after cleanup with restored c.
+    m.obj_val = dot(m.c, m.primal_sol) * (m.orig_sense == :Max ? -1 : +1)
 end
 
 status(m::ECOSMathProgModel) = m.solve_stat
@@ -391,8 +393,8 @@ function loadproblem!(m::ECOSMathProgModel, c, A, b, constr_cones, var_cones)
     m.ncones        = num_SOC_cones     # Num second-order cones
     m.conedims      = SOC_conedims      # Num contr. in each SOC
     m.nexp_cones    = num_exp_cones     # Num exponential cones
-    m.G             = ecos_G
-    m.A             = ecos_A
+    m.G             = ECOSMatrix(ecos_G)
+    m.A             = ECOSMatrix(ecos_A)
     m.c             = ecos_c
     m.orig_sense    = :Min
     m.h             = ecos_h