[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

@@ -89,6 +89,7 @@ extern "C"
 		restart_parameters_t restart_params;
 		double reflection_coefficient;
 		bool feasibility_polishing;
+		bool use_presolve;
 	} pdhg_parameters_t;
 
 	typedef struct

internal/presolve.h

@@ -0,0 +1,31 @@
+#ifndef PRESOLVE_H
+#define PRESOLVE_H
+
+#include "cupdlpx.h"
+#include "PSLP_API.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;
+} 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

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>
@@ -179,6 +180,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 +199,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 +234,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,167 @@
+#include "presolve.h"
+#include "cupdlpx.h"
+#include "PSLP_inf.h"
+#include "PSLP_sol.h"
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <float.h>
+#include <math.h>
+#include <time.h>
+
+void sanitize_infinity_for_pslp(double* arr, int n, double pslp_inf_val) {
+    for (int i = 0; i < n; ++i) {
+        if (isinf(arr[i]) || fabs(arr[i]) >= pslp_inf_val) {
+            arr[i] = (arr[i] > 0) ? pslp_inf_val : -pslp_inf_val;
+        }
+    }
+}
+
+#define CUPDLP_INF std::numeric_limits<double>::infinity()
+
+void restore_infinity_for_cupdlpx(double* arr, int n, double pslp_inf_val) {
+    for (int i = 0; i < n; ++i) {
+        if (arr[i] >= pslp_inf_val * 0.99) { 
+            arr[i] = INFINITY;
+        }
+        else if (arr[i] <= -pslp_inf_val * 0.99) {
+            arr[i] = -INFINITY;
+        }
+    }
+}
+
+lp_problem_t* convert_pslp_to_cupdlpx(PresolvedProblem *reduced_prob) {
+    matrix_desc_t desc;
+    memset(&desc, 0, sizeof(matrix_desc_t)); 
+
+    desc.m = reduced_prob->m;
+    desc.n = reduced_prob->n;
+    desc.fmt = matrix_csr;
+    desc.zero_tolerance = 0.0; 
+
+    desc.data.csr.nnz = reduced_prob->nnz;
+    desc.data.csr.row_ptr = reduced_prob->Ap;
+    desc.data.csr.col_ind = reduced_prob->Ai;
+    desc.data.csr.vals = reduced_prob->Ax;
+
+    double obj_const = 0.0; 
+
+    restore_infinity_for_cupdlpx(reduced_prob->lhs, reduced_prob->m, PSLP_INF);
+    restore_infinity_for_cupdlpx(reduced_prob->rhs, reduced_prob->m, PSLP_INF);
+    restore_infinity_for_cupdlpx(reduced_prob->lbs, reduced_prob->n, PSLP_INF);
+    restore_infinity_for_cupdlpx(reduced_prob->ubs, reduced_prob->n, PSLP_INF);
+
+    lp_problem_t* gpu_prob = create_lp_problem(
+        reduced_prob->c,
+        &desc,
+        reduced_prob->lhs, // constraint lower bound
+        reduced_prob->rhs, // constraint upper bound
+        reduced_prob->lbs, // variable lower bound
+        reduced_prob->ubs, // variable upper bound
+        &obj_const
+    );
+
+    return gpu_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. Sanitize input data
+    sanitize_infinity_for_pslp(original_prob->constraint_lower_bound, original_prob->num_constraints, PSLP_INF);
+    sanitize_infinity_for_pslp(original_prob->constraint_upper_bound, original_prob->num_constraints, PSLP_INF);
+    sanitize_infinity_for_pslp(original_prob->variable_lower_bound, original_prob->num_variables, PSLP_INF);
+    sanitize_infinity_for_pslp(original_prob->variable_upper_bound, original_prob->num_variables, PSLP_INF);
+
+    // 2. Init Settings
+    info->settings = default_settings();
+    if (params->verbose) {
+        info->settings->verbose = true;
+    }
+
+    // 3. 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 
+    );
+
+    // 4. 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->reduced_problem->objective_constant = original_prob->objective_constant;
+    }
+
+    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));
+        res->termination_reason = TERMINATION_REASON_PRIMAL_INFEASIBLE; // Simplified
+        res->cumulative_time_sec = info->presolve_time;
+        return res;
+    }
+
+    if (!reduced_result || !info->presolver) return NULL;
+
+    int n_red = info->presolver->reduced_prob->n;
+    double *z_dummy = (double*)calloc(n_red, sizeof(double)); 
+
+    postsolve(info->presolver, 
+              reduced_result->primal_solution, 
+              reduced_result->dual_solution, 
+              z_dummy, 
+              reduced_result->primal_objective_value);
+
+    free(z_dummy);
+
+    cupdlpx_result_t *final_result = (cupdlpx_result_t*)malloc(sizeof(cupdlpx_result_t));
+    *final_result = *reduced_result; 
+
+    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);
+}

src/solver.cu

@@ -19,6 +19,7 @@ limitations under the License.
 #include "preconditioner.h"
 #include "solver.h"
 #include "utils.h"
+#include "presolve.h"
 #include <cublas_v2.h>
 #include <cuda_runtime.h>
 #include <cusparse.h>
@@ -97,9 +98,30 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params,
                            const lp_problem_t *original_problem)
 {
     print_initial_info(params, original_problem);
-    rescale_info_t *rescale_info = rescale_problem(params, original_problem);
-    pdhg_solver_state_t *state =
-        initialize_solver_state(original_problem, rescale_info);
+
+    cupdlpx_presolve_info_t *presolve_info = NULL;
+    const lp_problem_t *problem_to_solve = original_problem;
+
+    if (params->use_presolve) {
+        if (params->verbose) printf("Running Presolve...\n");
+
+        presolve_info = pslp_presolve(original_problem, params);
+
+        if (!presolve_info) {
+             fprintf(stderr, "Presolve failed.\n");
+             return NULL;
+        }
+
+        if (presolve_info->problem_solved_during_presolve) {
+             cupdlpx_result_t *early_result = pslp_postsolve(presolve_info, NULL, original_problem);
+             cupdlpx_presolve_info_free(presolve_info);
+             return early_result;
+        }
+
+        problem_to_solve = presolve_info->reduced_problem;
+    }
+    rescale_info_t *rescale_info = rescale_problem(params, problem_to_solve);
+    pdhg_solver_state_t *state = initialize_solver_state(problem_to_solve, rescale_info);
 
     rescale_info_free(rescale_info);
     initialize_step_size_and_primal_weight(state, params);
@@ -164,9 +186,22 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params,
         feasibility_polish(params, state);
     }
 
-    cupdlpx_result_t *results = create_result_from_state(state);
+    cupdlpx_result_t *current_result = create_result_from_state(state);
+
+    cupdlpx_result_t *final_result = NULL;
+
+    if (params->use_presolve && presolve_info) {
+        final_result = pslp_postsolve(presolve_info, current_result, original_problem);
+
+        cupdlpx_result_free(current_result);
+        cupdlpx_presolve_info_free(presolve_info);
+    } 
+    else {
+        final_result = current_result;
+    }
+
     pdhg_solver_state_free(state);
-    return results;
+    return final_result;
 }
 
 static pdhg_solver_state_t *
@@ -977,6 +1012,8 @@ void set_default_parameters(pdhg_parameters_t *params)
     params->restart_params.k_i = 0.01;
     params->restart_params.k_d = 0.0;
     params->restart_params.i_smooth = 0.3;
+
+    params->use_presolve = false;
 }
 
 //Feasibility Polishing