[PRESOLVE] Avoid copying some rows/cols, fix issue unconstrained reduction zero with objective (#97)

* free memory, fix issue unconstrained reduction zero obj

* accelerate computation of presolved A and Q

* fix nnz errors, multipliers to dense vectors

* remove dropzeros!

* multipliers to dense vectors

* QuadraticModelPresolvedData

* fixed_vars_only

* up doc, rm setindex2

Author: geoffroyleconte

Project.toml

@@ -19,7 +19,7 @@ NLPModels = "0.18, 0.19"
 NLPModelsModifiers = "0.6"
 Requires = "0.3, 0.4, 0.5, 1.0"
 SolverCore = "0.2.2"
-SparseMatricesCOO = "0.1"
+SparseMatricesCOO = "0.1.1"
 julia = "^1.0.0"
 
 [extras]

src/presolve/empty_rows.jl

@@ -3,23 +3,21 @@ struct EmptyRow{T, S} <: PresolveOperation{T, S}
 end
 
 function empty_rows!(
+  qmp::QuadraticModelPresolveData{T, S},
   operations::Vector{PresolveOperation{T, S}},
-  lcon::S,
-  ucon::S,
-  ncon,
-  row_cnt,
-  kept_rows::Vector{Bool},
 ) where {T, S}
-  empty_row_pass = false
-  for i = 1:ncon
+  qmp.empty_row_pass = false
+  lcon, ucon, row_cnt, kept_rows = qmp.lcon, qmp.ucon, qmp.row_cnt, qmp.kept_rows
+  for i = 1:qmp.ncon
     (kept_rows[i] && (row_cnt[i] == 0)) || continue
-    empty_row_pass = true
+    qmp.empty_row_pass = true
     @assert (lcon[i] ≤ zero(T) ≤ ucon[i])
     row_cnt[i] = -1
     kept_rows[i] = false
     push!(operations, EmptyRow{T, S}(i))
   end
-  return empty_row_pass
 end
 
-postsolve!(sol::QMSolution, operation::EmptyRow) = nothing
+function postsolve!(sol::QMSolution, operation::EmptyRow, psd::PresolvedData)
+  psd.kept_rows[operation.i] = true
+end

src/presolve/free_rows.jl

@@ -3,22 +3,21 @@ struct FreeRow{T, S} <: PresolveOperation{T, S}
 end
 
 function free_rows!(
+  qmp::QuadraticModelPresolveData{T, S},
   operations::Vector{PresolveOperation{T, S}},
-  lcon::S,
-  ucon::S,
-  ncon,
-  row_cnt,
-  kept_rows::Vector{Bool},
 ) where {T, S}
-  free_row_pass = false
-  for i = 1:ncon
+  lcon, ucon = qmp.lcon, qmp.ucon
+  row_cnt, kept_rows = qmp.row_cnt, qmp.kept_rows
+  qmp.free_row_pass = false
+  for i = 1:qmp.ncon
     (kept_rows[i] && lcon[i] == -T(Inf) && ucon[i] == T(Inf)) || continue
-    free_row_pass = true
+    qmp.free_row_pass = true
     row_cnt[i] = -1
     kept_rows[i] = false
     push!(operations, FreeRow{T, S}(i))
   end
-  return free_row_pass
 end
 
-postsolve!(sol::QMSolution, operation::FreeRow) = nothing
+function postsolve!(sol::QMSolution, operation::FreeRow, psd::PresolvedData)
+  psd.kept_rows[operation.i] = true
+end

src/presolve/linear_singleton_columns.jl

@@ -2,92 +2,82 @@ struct FreeLinearSingletonColumn{T, S} <: PresolveOperation{T, S}
   i::Int
   j::Int
   aij::T
-  arowi::Row{T}
   yi::T
   conival::T
 end
 
 function free_linear_singleton_columns!(
+  qmp::QuadraticModelPresolveData{T, S},
   operations::Vector{PresolveOperation{T, S}},
-  hcols::Vector{Col{T}},
-  arows::Vector{Row{T}},
-  acols::Vector{Col{T}},
-  c::AbstractVector{T},
-  c0::T,
-  lcon::AbstractVector{T},
-  ucon::AbstractVector{T},
-  lvar::AbstractVector{T},
-  uvar::AbstractVector{T},
-  nvar,
-  row_cnt,
-  col_cnt,
-  kept_rows,
-  kept_cols,
 ) where {T, S}
-  free_lsc_pass = false
+  qmp.free_lsc_pass = false
   c0_offset = zero(T)
+  hcols, arows, acols, c = qmp.hcols, qmp.arows, qmp.acols, qmp.c
+  lcon, ucon, lvar, uvar = qmp.lcon, qmp.ucon, qmp.lvar, qmp.uvar
+  row_cnt, col_cnt = qmp.row_cnt, qmp.col_cnt
+  kept_rows, kept_cols = qmp.kept_rows, qmp.kept_cols
 
-  for j = 1:nvar
+  for j = 1:qmp.nvar
     (kept_cols[j] && (col_cnt[j] == 1)) || continue
     # check infinity bounds and no hessian contribution
     if lvar[j] == -T(Inf) && uvar[j] == T(Inf) && isempty(hcols[j].nzind)
-      # find i the row index of the singleton column j, and Aij
-      Aij = T(Inf)
+      # find i the row index of the singleton column j, and aij
+      aij = T(Inf)
       colj = acols[j]
       k = 1
       i = colj.nzind[k]
       while !(kept_rows[i])
         k += 1
         i = colj.nzind[k]
       end
-      Aij = colj.nzval[k]
+      aij = colj.nzval[k]
 
-      yi = c[j] / Aij
+      yi = c[j] / aij
       nzcj = c[j] != zero(T)
       if yi < zero(T)
         lcon[i] = ucon[i]
       elseif yi > zero(T)
         ucon[i] = lcon[i]
       end
-      if abs(Aij) > sqrt(eps(T))
+      if abs(aij) > sqrt(eps(T))
         nzcj && (c0_offset += yi * ucon[i]) # update c0
-        nb_elem_i = row_cnt[i] - 1
         rowi = arows[i] # i-th row
         # new row to store for postsolve:
-        rowi2 = Row(zeros(Int, nb_elem_i), zeros(T, nb_elem_i))
-        c_i = 1
-        for k = 1:length(rowi.nzind)
-          j2 = rowi.nzind[k]
+        for (l, ail) in zip(rowi.nzind, rowi.nzval)
           # add all col elements to rowi2 except the col j2 == j
-          if kept_cols[j2] && j2 != j
-            rowi2.nzind[c_i] = j2
-            Aij2 = rowi.nzval[k]
-            rowi2.nzval[c_i] = Aij2
-            nzcj && (c[j2] -= yi * Aij2) # update c if c[j] != 0
-            col_cnt[j2] -= 1
-            c_i += 1
-          end
+          (kept_cols[l] && l != j) || continue
+          nzcj && (c[l] -= yi * ail) # update c if c[j] != 0
+          col_cnt[l] -= 1
         end
         conival = nzcj ? ucon[i] : (lcon[i] + ucon[i]) / 2 # constant for postsolve
-        push!(operations, FreeLinearSingletonColumn{T, S}(i, j, Aij, rowi2, yi, conival))
+        push!(operations, FreeLinearSingletonColumn{T, S}(i, j, aij, yi, conival))
         kept_cols[j] = false
         col_cnt[j] = -1
         kept_rows[i] = false
         row_cnt[i] = -1
-        free_lsc_pass = true
+        qmp.free_lsc_pass = true
       end
     end
   end
-  return free_lsc_pass, c0 + c0_offset
+  qmp.c0 += c0_offset
 end
 
-function postsolve!(sol::QMSolution{T, S}, operation::FreeLinearSingletonColumn{T, S}) where {T, S}
+function postsolve!(
+  sol::QMSolution,
+  operation::FreeLinearSingletonColumn{T, S},
+  psd::PresolvedData{T, S},
+) where {T, S}
   x = sol.x
-  j = operation.j
-  # x[j] = (coival - Σₖ Aik x[k]) / Aij , where k ≂̸ j
+  kept_rows, kept_cols = psd.kept_rows, psd.kept_cols
+  i, j = operation.i, operation.j
+  arowi = psd.arows[i]
+  kept_rows[i] = true
+  kept_cols[j] = true
+  # x[j] = (coival - Σₖ Aik x[k]) / aij , where k ≂̸ j
   x[j] = operation.conival
-  for (i, Aij) in zip(operation.arowi.nzind, operation.arowi.nzval)
-    x[j] -= Aij * x[i]
+  for (l, ail) in zip(arowi.nzind, arowi.nzval)
+    (l != j && kept_cols[l]) || continue 
+    x[j] -= ail * x[l]
   end
   x[j] /= operation.aij
   sol.s_l[j] = zero(T)

src/presolve/postsolve_utils.jl

@@ -41,3 +41,28 @@ function restore_ilow_iupp!(ilow, iupp, kept_cols)
     end
   end
 end
+
+function restore_s!(
+  s_l::AbstractVector{T},
+  s_u::AbstractVector{T},
+  s_l_in::SparseVector{T, Int},
+  s_u_in::SparseVector{T, Int},
+  kept_cols::Vector{Bool},
+) where {T}
+  ilow, iupp = copy(s_l_in.nzind), copy(s_u_in.nzind)
+  restore_ilow_iupp!(ilow, iupp, kept_cols)
+  s_l[ilow] .= s_l_in.nzval
+  s_u[iupp] .= s_u_in.nzval
+end
+
+function restore_s!(
+  s_l::S,
+  s_u::S,
+  s_l_in::S,
+  s_u_in::S,
+  kept_cols::Vector{Bool},
+) where {S <: DenseVector}
+  nvar = length(s_l)
+  restore_x!(kept_cols, s_l_in, s_l, nvar)
+  restore_x!(kept_cols, s_u_in, s_u, nvar)
+end