[WIP] Integrate PSLP Presolver

CMakeLists.txt

@@ -72,6 +72,17 @@ if (CUPDLPX_BUILD_PYTHON)
     find_package(Python3 COMPONENTS Interpreter REQUIRED) # For versioning script and pybind11
 endif()
 
+include(FetchContent)
+
+FetchContent_Declare(
+  pslp
+  GIT_REPOSITORY https://github.com/dance858/PSLP.git
+  GIT_TAG        main
+)
+
+FetchContent_MakeAvailable(pslp)
+include_directories(${pslp_SOURCE_DIR}/PSLP)
+
 # -----------------------------------------------------------------------------
 # SOURCE DISCOVERY & TARGET DEFINITION
 # -----------------------------------------------------------------------------
@@ -97,6 +108,7 @@ set(CORE_LINK_LIBS PUBLIC
   CUDA::cublas
   CUDA::cusparse
   ZLIB::ZLIB
+  PSLP
 )
 
 # -----------------------------------------------------------------------------

include/cupdlpx_types.h

@@ -16,6 +16,7 @@ limitations under the License.
 
 #pragma once
 
+#include "PSLP_stats.h"
 #include <stdbool.h>
 #include <stdlib.h>
 
@@ -89,6 +90,7 @@ extern "C"
 		restart_parameters_t restart_params;
 		double reflection_coefficient;
 		bool feasibility_polishing;
+		bool use_presolve;
 	} pdhg_parameters_t;
 
 	typedef struct
@@ -98,10 +100,13 @@ extern "C"
 
 		double *primal_solution;
 		double *dual_solution;
+		double *reduced_cost;
 
 		int total_count;
 		double rescaling_time_sec;
 		double cumulative_time_sec;
+		double presolve_time;
+		double presolve_setup_time;
 
 		double absolute_primal_residual;
 		double relative_primal_residual;
@@ -118,6 +123,7 @@ extern "C"
 		termination_reason_t termination_reason;
 		double feasibility_polishing_time;
 		int feasibility_iteration;
+		PresolveStats presolve_stats;
 	} cupdlpx_result_t;
 
 	// matrix formats

internal/presolve.h

@@ -0,0 +1,33 @@
+#ifndef PRESOLVE_H
+#define PRESOLVE_H
+
+#include "cupdlpx.h"
+#include "PSLP_API.h"
+#include "PSLP_stats.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct {
+    Presolver *presolver;
+    Settings *settings;
+    lp_problem_t *reduced_problem; 
+    bool problem_solved_during_presolve; 
+    double presolve_time;
+    double presolve_setup_time;
+} cupdlpx_presolve_info_t;
+
+cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const pdhg_parameters_t *params);
+
+cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 
+                                 cupdlpx_result_t *reduced_result, 
+                                 const lp_problem_t *original_prob);
+
+void cupdlpx_presolve_info_free(cupdlpx_presolve_info_t *info);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif // PRESOLVE_H

internal/utils.h

@@ -99,10 +99,14 @@ extern "C"
 
     void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem);
 
-    void pdhg_final_log(
-        const pdhg_solver_state_t *solver_state,
-        bool verbose,
-        termination_reason_t termination_reason);
+    // void pdhg_final_log(
+    //     const pdhg_solver_state_t *solver_state,
+    //     bool verbose,
+    //     termination_reason_t termination_reason);
+    void pdhg_final_log(const pdhg_solver_state_t *state, 
+                    const PresolveStats *stats,
+                    bool verbose,
+                    termination_reason_t reason);
 
     void display_iteration_stats(const pdhg_solver_state_t *solver_state, bool verbose);
 

src/cli.c

@@ -18,6 +18,7 @@ limitations under the License.
 #include "mps_parser.h"
 #include "solver.h"
 #include "utils.h"
+#include "presolve.h"
 #include <getopt.h>
 #include <libgen.h>
 #include <stdbool.h>
@@ -139,6 +140,33 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir,
         fprintf(outfile, "Feasibility Polishing Time (sec): %e\n", result->feasibility_polishing_time);
         fprintf(outfile, "Feasibility Polishing Iteration Count: %d\n", result->feasibility_iteration);
     }
+    if (result->presolve_time > 0.0) {
+        fprintf(outfile, "Presolve Time (sec): %e\n", result->presolve_time);
+        fprintf(outfile, "Presolve Setup Time (sec): %e\n", result->presolve_setup_time);
+    }
+    if (result->presolve_stats.n_cols_original > 0) {
+        fprintf(outfile, "Original Rows: %d\n", result->presolve_stats.n_rows_original);
+        fprintf(outfile, "Original Cols: %d\n", result->presolve_stats.n_cols_original);
+        fprintf(outfile, "Original NNZ: %d\n", result->presolve_stats.nnz_original);
+        fprintf(outfile, "Reduced Rows: %d\n", result->presolve_stats.n_rows_reduced);
+        fprintf(outfile, "Reduced Cols: %d\n", result->presolve_stats.n_cols_reduced);
+        fprintf(outfile, "Reduced NNZ: %d\n", result->presolve_stats.nnz_reduced);
+
+        fprintf(outfile, "NNZ Removed Trivial: %d\n", result->presolve_stats.nnz_removed_trivial);
+        fprintf(outfile, "NNZ Removed Fast: %d\n", result->presolve_stats.nnz_removed_fast);
+        fprintf(outfile, "NNZ Removed Primal Propagation: %d\n", result->presolve_stats.nnz_removed_primal_propagation);
+        fprintf(outfile, "NNZ Removed Parallel Rows: %d\n", result->presolve_stats.nnz_removed_parallel_rows);
+        fprintf(outfile, "NNZ Removed Parallel Cols: %d\n", result->presolve_stats.nnz_removed_parallel_cols);
+
+        fprintf(outfile, "Presolve Total Time (sec): %e\n", result->presolve_stats.presolve_total_time);
+        fprintf(outfile, "Presolve Time Init (sec): %e\n", result->presolve_stats.ps_time_init);
+        fprintf(outfile, "Presolve Time Fast (sec): %e\n", result->presolve_stats.ps_time_fast);
+        fprintf(outfile, "Presolve Time Medium (sec): %e\n", result->presolve_stats.ps_time_medium);
+        fprintf(outfile, "Presolve Time Primal Proppagation (sec): %e\n", result->presolve_stats.ps_time_primal_propagation);
+        fprintf(outfile, "Presolve Time Parallel Rows (sec): %e\n", result->presolve_stats.ps_time_parallel_rows);
+        fprintf(outfile, "Presolve Time Parallel Cols (sec): %e\n", result->presolve_stats.ps_time_parallel_cols);
+        fprintf(outfile, "Postsolve Time (sec): %e\n", result->presolve_stats.ps_time_post_solve);
+    }
     fclose(outfile);
     free(file_path);
 }
@@ -179,6 +207,8 @@ void print_usage(const char *prog_name)
                     "polish tolerance (default: 1e-6).\n");
     fprintf(stderr, "  -f  --feasibility_polishing         Enable feasibility " 
                     "use feasibility polishing (default: false).\n");
+    fprintf(stderr, "  -p, --presolve "
+                    "enable presolving (default: false).\n");
 }
 
 int main(int argc, char *argv[])
@@ -196,10 +226,11 @@ int main(int argc, char *argv[])
         {"eps_infeas_detect", required_argument, 0, 1005},
         {"eps_feas_polish", required_argument, 0, 1006},
         {"feasibility_polishing", no_argument, 0, 'f'},
+        {"presolve", no_argument, 0, 'p'},
         {0, 0, 0, 0}};
 
     int opt;
-    while ((opt = getopt_long(argc, argv, "hvf", long_options, NULL)) != -1)
+    while ((opt = getopt_long(argc, argv, "hvfp", long_options, NULL)) != -1)
     {
         switch (opt)
         {
@@ -230,6 +261,9 @@ int main(int argc, char *argv[])
         case 'f':                  // --feasibility_polishing
             params.feasibility_polishing = true;
             break;
+        case 'p':                  // --presolve
+            params.use_presolve = true;
+            break;
         case '?': // Unknown option
             return 1;
         }

src/presolve.c

@@ -0,0 +1,148 @@
+#include "presolve.h"
+#include "cupdlpx.h"
+#include "PSLP_sol.h"
+#include "utils.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <time.h>
+
+lp_problem_t* convert_pslp_to_cupdlpx(PresolvedProblem *reduced_prob) {
+
+    lp_problem_t *cupdlpx_prob = (lp_problem_t *)safe_malloc(sizeof(lp_problem_t));
+    // TODO: handle warmstart here
+    cupdlpx_prob->primal_start = NULL;
+    cupdlpx_prob->dual_start = NULL;
+
+    cupdlpx_prob->objective_constant = reduced_prob->obj_offset;
+    cupdlpx_prob->objective_vector = reduced_prob->c;
+
+    cupdlpx_prob->constraint_lower_bound = reduced_prob->lhs;
+    cupdlpx_prob->constraint_upper_bound = reduced_prob->rhs;
+    cupdlpx_prob->variable_lower_bound = reduced_prob->lbs;
+    cupdlpx_prob->variable_upper_bound = reduced_prob->ubs;
+
+    cupdlpx_prob->constraint_matrix_num_nonzeros = reduced_prob->nnz;
+    cupdlpx_prob->constraint_matrix_row_pointers = reduced_prob->Ap;
+    cupdlpx_prob->constraint_matrix_col_indices = reduced_prob->Ai;
+    cupdlpx_prob->constraint_matrix_values = reduced_prob->Ax;
+
+    cupdlpx_prob->num_variables = reduced_prob->n;
+    cupdlpx_prob->num_constraints = reduced_prob->m;
+
+    return cupdlpx_prob;
+}
+
+cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const pdhg_parameters_t *params) {
+    clock_t start_time = clock();
+
+    cupdlpx_presolve_info_t *info = (cupdlpx_presolve_info_t*)calloc(1, sizeof(cupdlpx_presolve_info_t));
+    if (!info) return NULL;
+
+    // 1. Init Settings
+    info->settings = default_settings();
+    if (params->verbose) {
+        info->settings->verbose = true;
+    }
+    // info->settings->relax_bounds = false;
+
+    // 2. Init Presolver
+    info->presolver = new_presolver(
+        original_prob->constraint_matrix_values,
+        original_prob->constraint_matrix_col_indices,
+        original_prob->constraint_matrix_row_pointers,
+        original_prob->num_constraints,
+        original_prob->num_variables,
+        original_prob->constraint_matrix_num_nonzeros,
+        original_prob->constraint_lower_bound,
+        original_prob->constraint_upper_bound,
+        original_prob->variable_lower_bound,
+        original_prob->variable_upper_bound,
+        original_prob->objective_vector,
+        info->settings,
+        true 
+    );
+    info->presolve_setup_time = (double)(clock() - start_time) / CLOCKS_PER_SEC;
+
+    // 3. Run Presolve
+    PresolveStatus status = run_presolver(info->presolver);
+
+    if (status & INFEASIBLE || status & UNBOUNDED) {
+        info->problem_solved_during_presolve = true;
+        info->reduced_problem = NULL;
+        if (params->verbose) printf("Problem solved by presolver (Infeasible/Unbounded).\n");
+    } else {
+        info->problem_solved_during_presolve = false;
+        if (params->verbose) {
+            printf("Presolve finished. Reduced rows: %d, cols: %d\n", 
+                   info->presolver->reduced_prob->m, info->presolver->reduced_prob->n);
+        }
+        info->reduced_problem = convert_pslp_to_cupdlpx(info->presolver->reduced_prob);
+    }
+
+    info->presolve_time = (double)(clock() - start_time) / CLOCKS_PER_SEC;
+    return info;
+}
+
+cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 
+                                 cupdlpx_result_t *reduced_result, 
+                                 const lp_problem_t *original_prob) {
+
+    if (info->problem_solved_during_presolve) {
+        cupdlpx_result_t *res = (cupdlpx_result_t*)calloc(1, sizeof(cupdlpx_result_t));
+
+        if (info->presolver && info->presolver->stats) {
+            res->presolve_stats = *(info->presolver->stats);
+        }
+
+        res->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE;
+        res->cumulative_time_sec = info->presolve_time;
+        res->num_variables = original_prob->num_variables;
+        res->num_constraints = original_prob->num_constraints;
+
+        if (res->num_variables > 0) {
+            res->primal_solution = (double*)calloc(res->num_variables, sizeof(double));
+            res->reduced_cost = (double*)calloc(res->num_variables, sizeof(double)); 
+        }
+        if (res->num_constraints > 0) {
+            res->dual_solution = (double*)calloc(res->num_constraints, sizeof(double));
+        }
+        return res;
+    }
+
+    if (!reduced_result || !info->presolver) return NULL;
+
+    postsolve(info->presolver, 
+              reduced_result->primal_solution, 
+              reduced_result->dual_solution, 
+              reduced_result->reduced_cost,
+              reduced_result->primal_objective_value);
+
+    cupdlpx_result_t *final_result = (cupdlpx_result_t*)malloc(sizeof(cupdlpx_result_t));
+
+    *final_result = *reduced_result; 
+
+    if (info->presolver->stats != NULL) {
+        final_result->presolve_stats = *(info->presolver->stats);
+    }
+
+    final_result->primal_solution = (double*)malloc(original_prob->num_variables * sizeof(double));
+    final_result->dual_solution = (double*)malloc(original_prob->num_constraints * sizeof(double));
+
+    memcpy(final_result->primal_solution, info->presolver->sol->x, original_prob->num_variables * sizeof(double));
+    memcpy(final_result->dual_solution, info->presolver->sol->y, original_prob->num_constraints * sizeof(double));
+    final_result->primal_objective_value = info->presolver->sol->obj;
+    final_result->cumulative_time_sec += info->presolve_time;
+
+    return final_result;
+}
+
+void cupdlpx_presolve_info_free(cupdlpx_presolve_info_t *info) {
+    if (!info) return;
+    // if (info->reduced_problem) lp_problem_free(info->reduced_problem);
+    if (info->presolver) free_presolver(info->presolver);
+    if (info->settings) free_settings(info->settings);
+    free(info);
+}