Release large PSLP presolve buffers early

CMakeLists.txt

@@ -3,7 +3,7 @@ project(PSLP LANGUAGES C CXX) # The library is written in C, but we also have so
 
 set(PSLP_VERSION_MAJOR 0)
 set(PSLP_VERSION_MINOR 0)
-set(PSLP_VERSION_PATCH 10)
+set(PSLP_VERSION_PATCH 11)
 set(PSLP_VERSION "${PSLP_VERSION_MAJOR}.${PSLP_VERSION_MINOR}.${PSLP_VERSION_PATCH}")
 add_compile_definitions(PSLP_VERSION="${PSLP_VERSION}")
 

include/PSLP/PSLP_API.h

@@ -124,6 +124,14 @@ extern "C"
        presolver contains the presolved problem. */
     PresolveStatus run_presolver(Presolver *presolver);
 
+    /* Free the presolved problem data owned by 'presolver' after the caller has
+       copied all arrays in 'reduced_prob'. This does not free 'presolver'; the
+       caller must still call 'free_presolver' after postsolve. Normal postsolve
+       and dual-infeasibility-ray postsolve remain available. Primal-
+       infeasibility-ray postsolve is no longer available because it requires
+       the reduced constraint matrix. */
+    void free_presolver_reduced_problem(Presolver *presolver);
+
     /* Postsolve the problem given the primal-dual solution (x, y, z) of the
        reduced problem. The function populates presolver->sol, so if you're
        looking for the solution to the original problem, you want to look there.

src/core/Presolver.c

@@ -736,6 +736,12 @@ PresolveStatus run_presolver(Presolver *presolver)
     stats->nnz_reduced = A->nnz;
     DEBUG(run_debugger_stats_consistency_check(stats));
     populate_presolved_problem(presolver);
+
+    // The transpose is only needed during presolve. Postsolve uses the recorded
+    // reductions and row/column mappings.
+    free_matrix(prob->constraints->AT);
+    prob->constraints->AT = NULL;
+
     clock_gettime(CLOCK_MONOTONIC, &outer_timer.end);
     stats->time_presolve = GET_ELAPSED_SECONDS(outer_timer);
 
@@ -747,6 +753,50 @@ PresolveStatus run_presolver(Presolver *presolver)
     return final_status(stats);
 }
 
+void free_presolver_reduced_problem(Presolver *presolver)
+{
+
+    // 'Ap' is non-null when 'run_presolver' has populated the reduced problem.
+    // Bailing out here keeps this function safe to call before a successful
+    // presolve.
+    if (!presolver || !presolver->reduced_prob || !presolver->reduced_prob->Ap)
+    {
+        return;
+    }
+
+    PresolvedProblem *reduced_prob = presolver->reduced_prob;
+    Problem *prob = presolver->prob;
+    if (prob)
+    {
+        Constraints *constraints = prob->constraints;
+        if (constraints)
+        {
+            free_matrix(constraints->A);
+            constraints->A = NULL;
+            free_matrix(constraints->AT);
+            constraints->AT = NULL;
+            PS_FREE(constraints->lhs);
+            PS_FREE(constraints->rhs);
+            PS_FREE(constraints->bounds);
+            PS_FREE(constraints->row_tags);
+            PS_FREE(constraints->col_tags);
+        }
+        if (prob->obj)
+        {
+            PS_FREE(prob->obj->c);
+        }
+    }
+
+    reduced_prob->Ax = NULL;
+    reduced_prob->Ai = NULL;
+    reduced_prob->lhs = NULL;
+    reduced_prob->rhs = NULL;
+    reduced_prob->c = NULL;
+    PS_FREE(reduced_prob->Ap);
+    PS_FREE(reduced_prob->lbs);
+    PS_FREE(reduced_prob->ubs);
+}
+
 void postsolve(Presolver *presolver, const double *x, const double *y,
                const double *z)
 {
@@ -798,6 +848,14 @@ void postsolve_primal_infeas_ray(Presolver *presolver, const double *y,
     size_t len_z = MAX((size_t) 1, stats->n_cols_reduced);
     double *z;
     assert(reduced_prob != NULL);
+    if (!reduced_prob || !reduced_prob->Ax)
+    {
+        // The reduced constraint matrix was released by
+        // 'free_presolver_reduced_problem', so the ray cannot be recovered.
+        fprintf(stderr, "PSLP warning: postsolve_primal_infeas_ray called after "
+                        "free_presolver_reduced_problem! Don't do this! \n");
+        return;
+    }
     z = (double *) ps_malloc(len_z, sizeof(double));
     assert(z != NULL);
 

tests/test_postsolve.h

@@ -46,12 +46,26 @@ static char *test_0_postsolve()
     Constraints *constraints = prob->constraints;
     Matrix *A = constraints->A;
 
-    run_presolver(presolver);
+    // calling before presolve must be a no-op that keeps the original problem
+    free_presolver_reduced_problem(presolver);
+    mu_assert("free before presolve must not touch the problem",
+              constraints->A != NULL && constraints->AT != NULL);
 
+    run_presolver(presolver);
+    mu_assert("transpose should be released after presolve",
+              constraints->AT == NULL);
     Mapping *maps = prob->constraints->state->work->mappings;
     int *rows_map = maps->rows;
     int *cols_map = maps->cols;
 
+    // should be able to call it multiple times
+    free_presolver_reduced_problem(presolver);
+    free_presolver_reduced_problem(presolver);
+    mu_assert("reduced matrix should be released", constraints->A == NULL);
+    mu_assert("reduced problem buffers should be released",
+              presolver->reduced_prob->Ax == NULL &&
+                  presolver->reduced_prob->Ap == NULL);
+
     // construct optimal primal solution to reduced problem (computed offline)
     double x[] = {10., -10., -10., 2.71428571, -10., -6.85714286, -10.};
     double y[] = {-2.57142857, 0.14285714};

tests/test_ray_postsolve.h

@@ -242,6 +242,7 @@ static char *test_0_dual_ray_postsolve()
 
     PresolveStatus status = run_presolver(presolver);
     mu_assert("dual ray postsolve presolve status", status == UNCHANGED);
+    free_presolver_reduced_problem(presolver);
 
     double x[] = {0.0, 0.0};
     double x_orig[] = {0.0, 0.0};