empty rows

Author: geoffroyleconte

src/presolve/empty_rows.jl

@@ -0,0 +1,49 @@
+function row_cnt!(Arows, row_cnt::Vector{Int})
+  for k=1:length(Arows)
+    i = Arows[k]
+    row_cnt[i] += 1
+  end
+end
+
+function removed_rows(row_cnt::Vector{Int})
+  return findall(isequal(0), row_cnt)
+end
+
+function empty_rows!(Arows, lcon::Vector{T}, ucon::Vector{T}, ncon, row_cnt::Vector{Int}, 
+                     rows_rm::Vector{Int}, Arows_sortperm::Vector{Int}) where {T}
+  new_ncon = 0
+  for i=1:ncon
+    if row_cnt[i] == 0
+      @assert lcon[i] ≤ zero(T) ≤ ucon[i]
+    else
+      new_ncon += 1
+      lcon[new_ncon] = lcon[i]
+      ucon[new_ncon] = ucon[i]
+    end
+  end
+  resize!(lcon, new_ncon)
+  resize!(ucon, new_ncon)
+
+  Arows_s = @views Arows[Arows_sortperm]
+  c_rm = 1
+  nrm = length(rows_rm)
+  for k=1:length(Arows)
+    while c_rm ≤ nrm && Arows_s[k] ≥ rows_rm[c_rm]
+      c_rm += 1 
+    end
+    Arows_s[k] -= c_rm - 1
+  end
+  return new_ncon
+end
+
+function restore_y!(y::Vector{T}, yout::Vector{T}, row_cnt, ncon) where {T}
+  c_y = 0
+  for i = 1:ncon
+    if row_cnt[i] == 0
+      yout[i] = zero(T)
+    else
+      c_y += 1
+      yout[i] = y[c_y]
+    end
+  end
+end

src/presolve/presolve.jl

@@ -1,10 +1,17 @@
 include("remove_ifix.jl")
+include("empty_rows.jl")
+
+mutable struct PresolvedData{T, S}
+  xrm::S
+  row_cnt::Vector{Int}
+  nconps::Int
+end
 
 mutable struct PresolvedQuadraticModel{T, S, M1, M2} <: AbstractQuadraticModel{T, S}
   meta::NLPModelMeta{T, S}
   counters::Counters
   data::QPData{T, S, M1, M2}
-  xrm::S
+  psd::PresolvedData{T, S}
 end
 
 """
@@ -64,6 +71,13 @@ function presolve(
     xrm = S(undef, 0)
   end
 
+  # remove constraints
+  row_cnt = zeros(Int, ncon)
+  row_cnt!(psdata.A.rows, row_cnt)
+  rows_rm = removed_rows(row_cnt)
+  Arows_sortperm = sortperm(psdata.A.rows)
+  nconps = empty_rows!(psdata.A.rows, lcon, ucon, ncon, row_cnt, rows_rm, Arows_sortperm)
+
   # form meta
   nnzh = length(psdata.H.vals)
   if !(nnzh == length(psdata.H.rows) == length(psdata.H.cols))
@@ -86,7 +100,7 @@ function presolve(
       qm,
       solution = xrm,
       objective = obj(qm, xrm),
-      multipliers = zeros(T, qm.meta.nvar),
+      multipliers = zeros(T, nconps),
       multipliers_L = s_l,
       multipliers_U = s_u,
       iter = 0,
@@ -98,17 +112,18 @@ function presolve(
       nvarps,
       lvar = lvar,
       uvar = uvar,
-      ncon = ncon,
+      ncon = nconps,
       lcon = lcon,
       ucon = ucon,
       nnzj = nnzj,
       nnzh = nnzh,
-      lin = 1:ncon,
+      lin = 1:nconps,
       islp = (nnzh == 0);
       minimize = qm.meta.minimize,
       kwargs...,
     )
-    ps = PresolvedQuadraticModel(psmeta, Counters(), psdata, xrm)
+    psd = PresolvedData{T, S}(xrm, row_cnt, nconps)
+    ps = PresolvedQuadraticModel(psmeta, Counters(), psdata, psd)
     return GenericExecutionStats(
       :unknown,
       ps,
@@ -131,10 +146,15 @@ function postsolve!(
   psqm::PresolvedQuadraticModel{T, S},
   x_in::S,
   x_out::S,
+  y_in::S,
+  y_out::S,
 ) where {T, S}
   if length(qm.meta.ifix) > 0
-    restore_ifix!(qm.meta.ifix, psqm.xrm, x_in, x_out)
+    restore_ifix!(qm.meta.ifix, psqm.psd.xrm, x_in, x_out)
   else
     x_out .= @views x_in[1:(qm.meta.nvar)]
   end
+  ncon = length(y_out)
+  restore_y!(y_in, y_out, psqm.psd.row_cnt, ncon)
+
 end

src/qpmodel.jl

@@ -94,6 +94,8 @@ function QuadraticModel(
   sortcols::Bool = false,
   kwargs...,
 ) where {T, S}
+  @assert all(lvar .≤ uvar)
+  @assert all(lcon .≤ ucon)
   nnzh = length(Hvals)
   if !(nnzh == length(Hrows) == length(Hcols))
     error("The length of Hrows, Hcols and Hvals must be the same")
@@ -155,6 +157,8 @@ function QuadraticModel(
   c0::T = zero(T),
   kwargs...,
 ) where {T, S}
+  @assert all(lvar .≤ uvar)
+  @assert all(lcon .≤ ucon)
   ncon, nvar = size(A)
   if typeof(H) <: AbstractLinearOperator # convert A to a LinOp if A is a Matrix?
     nnzh = 0

test/test_presolve.jl

@@ -48,13 +48,15 @@
   @test norm(Aps - sparse(Aps_true)) ≤ sqrt(eps(T))
   @test psqp.meta.lvar == lvarps_true
   @test psqp.meta.uvar == uvarps_true
-  @test psqp.xrm == [2.0]
+  @test psqp.psd.xrm == [2.0]
   @test psqp.meta.ifix == Int[]
   @test psqp.meta.nvar == 2
 
   x_in = [4.0; 7.0]
   x_out = zeros(3)
-  postsolve!(qp, psqp, x_in, x_out)
+  y_in = [2.0; 2.0]
+  y_out = zeros(2)
+  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
   @test x_out == [4.0; 2.0; 7.0]
 
   # test that solves the problem
@@ -68,3 +70,57 @@
   stats_ps2 = presolve(qp2)
   @test stats_ps2.status == :first_order
 end
+
+@testset "presolve empty rows" begin
+  H = [
+    6.0 2.0 1.0
+    2.0 5.0 2.0
+    1.0 2.0 4.0
+  ]
+  c = [-8.0; -3; -3]
+  A = [
+    1.0 0.0 1.0
+    0.0 0.0 0.0
+    0.0 0.0 0.0
+    3.2 0.0 0.0
+    0.0 0.0 0.0
+    0.0 2.0 1.0
+  ]
+  b = [0.0; 0.0; 0.0; 4.0; 0.0; 3]
+  l = [0.0; 0.0; 0]
+  u = [Inf; 2.0; Inf]
+  T = eltype(c)
+  qp = QuadraticModel(
+    c,
+    SparseMatrixCOO(tril(H)),
+    A = SparseMatrixCOO(A),
+    lcon = b,
+    ucon = b,
+    lvar = l,
+    uvar = u,
+    c0 = 0.0,
+    name = "QM1",
+  )
+
+  statsps = presolve(qp)
+  psqp = statsps.solver_specific[:presolvedQM]
+
+  Aps_true = [
+    1.0 0.0 1.0
+    3.2 0.0 0.0
+    0.0 2.0 1.0
+  ]
+  bps_true = [0.0; 4.0; 3.0]
+
+  Aps = sparse(psqp.data.A.rows, psqp.data.A.cols, psqp.data.A.vals, psqp.meta.ncon, psqp.meta.nvar)
+  @test Aps == sparse(Aps_true)
+  @test psqp.meta.lcon == psqp.meta.ucon == bps_true
+  @test psqp.meta.ncon == 3
+
+  x_in = [4.0; 7.0; 4.0]
+  x_out = zeros(3)
+  y_in = [2.0; 2.0; 4.0]
+  y_out = zeros(6)
+  postsolve!(qp, psqp, x_in, x_out, y_in, y_out)
+  @test y_out == [2.0; 0.0; 0.0; 2.0; 0.0; 4.0]
+end