[PRESOLVE] Free linear singleton columns, row/col storage of A and H (#87)

* clean presolve

* up postsolve

* remove useless storage vectors

* store rows and cols

* free linear singleton columns

Author: geoffroyleconte

src/presolve/empty_rows.jl

@@ -1,51 +1,25 @@
-"""
-    new_ncon = empty_rows!(Arows, lcon, ucon, ncon, row_cnt, empty_rows)
-
-Removes the empty rows of A, and the corresponding elements in lcon and ucon that are in `empty_rows`.
-`row_cnt` is a vector of the number of elements per row.
+struct EmptyRow{T, S} <: PresolveOperation{T, S}
+  i::Int # idx of empty row
+end
 
-Returns the new number of constraints `new_ncon` and updates in-place `Arows`, `lcon`, `ucon`.
-"""
 function empty_rows!(
-  Arows,
-  lcon::Vector{T},
-  ucon::Vector{T},
+  operations::Vector{PresolveOperation{T, S}},
+  lcon::S,
+  ucon::S,
   ncon,
-  row_cnt::Vector{Int},
-  empty_rows::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
+  row_cnt,
+  kept_rows::Vector{Bool},
+) where {T, S}
+  empty_row_pass = false
+  for i=1:ncon
+    (kept_rows[i] && (row_cnt[i] == 0)) || continue
+    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
-  resize!(lcon, new_ncon)
-  resize!(ucon, new_ncon)
-
-  # assume Acols is sorted
-  Annz = length(Arows)
-  nempty = length(empty_rows)
-
-  for idxempty = 1:nempty
-    currentiempty = empty_rows[idxempty]
-    # index of the current singleton row that takes the number of 
-    # already removed variables into account:
-    newcurrentiempty = currentiempty - idxempty + 1
-
-    # remove singleton rows in A rows
-    Awritepos = 1
-    k = 1
-    while k <= Annz
-      Ai = Arows[k]
-      Arows[Awritepos] = (Ai < newcurrentiempty) ? Ai : Ai - 1
-      Awritepos += 1
-      k += 1
-    end
-  end
-  return new_ncon
+  return empty_row_pass
 end
+
+postsolve!(pt::OutputPoint, operation::EmptyRow) = nothing

src/presolve/linear_singleton_columns.jl

@@ -0,0 +1,96 @@
+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!(
+  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
+  c0_offset = zero(T)
+
+  for j=1: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)
+      colj = acols[j]
+      k = 1
+      i = colj.nzind[k]
+      while !(kept_rows[i])
+        k += 1
+        i = colj.nzind[k]
+      end
+      Aij = colj.nzval[k]
+
+      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))
+        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 in 1:length(rowi.nzind)
+          j2 = rowi.nzind[k]
+          # 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
+        end
+        conival = nzcj ? ucon[i] : (lcon[i] + ucon[i]) / 2 # constant for postsolve
+        push!(operations, FreeLinearSingletonColumn{T, S}(i, j, Aij, rowi2, yi, conival))
+        kept_cols[j] = false
+        col_cnt[j] = -1
+        kept_rows[i] = false
+        row_cnt[i] = -1
+        free_lsc_pass = true
+      end
+    end
+  end
+  return free_lsc_pass, c0 + c0_offset
+end
+
+function postsolve!(pt::OutputPoint{T, S}, operation::FreeLinearSingletonColumn{T, S}) where {T, S}
+  x = pt.x
+  j = operation.j
+  # 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]
+  end
+  x[j] /= operation.aij
+  pt.s_l[j] = zero(T)
+  pt.s_u[j] = zero(T)
+  pt.y[operation.i] = operation.yi
+end

src/presolve/postsolve_utils.jl

@@ -1,18 +1,15 @@
-function restore_ifix!(kept_cols, xps, x, xout)
+function restore_x!(kept_cols, xin::S, xout::S, nvar) where S
   # put x and xps inside xout according to kept_cols
-  nvar = length(xout)
   cx = 0
   for i = 1:nvar
     if kept_cols[i]
       cx += 1
-      xout[i] = x[cx]
-    else
-      xout[i] = xps[i]
+      xout[i] = xin[cx]
     end
   end
 end
 
-function restore_y!(y::Vector{T}, yout::Vector{T}, kept_rows::Vector{Bool}, ncon) where {T}
+function restore_y!(kept_rows::Vector{Bool}, y::Vector{T}, yout::Vector{T}, ncon) where {T}
   c_y = 0
   for i = 1:ncon
     if !kept_rows[i]