Integrate PSLP Presolver (#41)

* add PSLP presolver

* fix: correct the reduced cost in PSLP postsolve

* add presolve stats

* test PSLP relax_bounds false

* remove PSLP_INF

* add presolve_setup_time and presolve_time record, improve reduce prob memory

* set relax_bounds to True (default)

* fix Segmentation fault when presolve detected infeasible

* refactor: PSLP presolve

* add verbose check for log

* set use_presolve to true by default

* rename use_presolve to presolve

* refactor: add presolve in python interface

* docs: add presolve in readme

* feat: release v0.2.0

* fix: off-by-one iteration error and ensure consistent final metrics

Author: ZedongPeng

CMakeLists.txt

@@ -7,8 +7,8 @@ cmake_minimum_required(VERSION 3.20)
 project(cupdlpx LANGUAGES C CXX CUDA)
 
 set(CUPDLPX_VERSION_MAJOR 0)
-set(CUPDLPX_VERSION_MINOR 1)
-set(CUPDLPX_VERSION_PATCH 5)
+set(CUPDLPX_VERSION_MINOR 2)
+set(CUPDLPX_VERSION_PATCH 0)
 
 set(CUPDLPX_VERSION "${CUPDLPX_VERSION_MAJOR}.${CUPDLPX_VERSION_MINOR}.${CUPDLPX_VERSION_PATCH}")
 add_compile_definitions(CUPDLPX_VERSION="${CUPDLPX_VERSION}")
@@ -75,6 +75,20 @@ if (CUPDLPX_BUILD_PYTHON)
     find_package(Python3 COMPONENTS Interpreter REQUIRED) # For versioning script and pybind11
 endif()
 
+include(FetchContent)
+
+set(PSLP_VERSION_TAG "v0.0.2")
+
+FetchContent_Declare(
+  pslp
+  GIT_REPOSITORY https://github.com/dance858/PSLP.git
+  GIT_TAG        ${PSLP_VERSION_TAG}
+)
+
+FetchContent_MakeAvailable(pslp)
+include_directories(${pslp_SOURCE_DIR}/PSLP)
+add_compile_definitions(PSLP_VERSION=\"${PSLP_VERSION_TAG}\")
+
 # -----------------------------------------------------------------------------
 # SOURCE DISCOVERY & TARGET DEFINITION
 # -----------------------------------------------------------------------------
@@ -100,6 +114,7 @@ set(CORE_LINK_LIBS PUBLIC
   CUDA::cublas
   CUDA::cusparse
   ZLIB::ZLIB
+  PSLP
 )
 
 # -----------------------------------------------------------------------------

README.md

@@ -94,6 +94,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the
 | `--eval_freq` | `int` | Termination evaluation frequency | `200` |
 | `--sv_max_iter` | `int` | Max iterations for singular value estimation | `5000` |
 | `--sv_tol` | `float` | Tolerance for singular value estimation | `1e-4` |
+| `--no_presolve` | `flag` | Disable presolve | `enabled` |
 | `-f`,`--feasibility_polishing` |`flag` | Run the polishing loop | `false` |
 | `--eps_feas_polish` | `double` | Relative tolerance for polishing | `1e-6`  |
 

include/cupdlpx_types.h

@@ -16,6 +16,7 @@ limitations under the License.
 
 #pragma once
 
+#include "PSLP_stats.h"
 #include <stdbool.h>
 #include <stdlib.h>
 
@@ -30,6 +31,7 @@ extern "C"
 		TERMINATION_REASON_OPTIMAL,
 		TERMINATION_REASON_PRIMAL_INFEASIBLE,
 		TERMINATION_REASON_DUAL_INFEASIBLE,
+		TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED,
 		TERMINATION_REASON_TIME_LIMIT,
 		TERMINATION_REASON_ITERATION_LIMIT,
 		TERMINATION_REASON_FEAS_POLISH_SUCCESS
@@ -91,19 +93,28 @@ extern "C"
 		restart_parameters_t restart_params;
 		double reflection_coefficient;
 		bool feasibility_polishing;
+		bool presolve;
 	} pdhg_parameters_t;
 
 	typedef struct
 	{
 		int num_variables;
 		int num_constraints;
+		int num_nonzeros;
+
+		int num_reduced_variables;
+		int num_reduced_constraints;
+		int num_reduced_nonzeros;
 
 		double *primal_solution;
 		double *dual_solution;
+		double *reduced_cost;
 
 		int total_count;
 		double rescaling_time_sec;
 		double cumulative_time_sec;
+		double presolve_time;
+		int presolve_status;
 
 		double absolute_primal_residual;
 		double relative_primal_residual;

internal/presolve.h

@@ -0,0 +1,39 @@
+#ifndef PRESOLVE_H
+#define PRESOLVE_H
+
+#include "PSLP_API.h"
+#include "PSLP_stats.h"
+#include "cupdlpx.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;
+        char presolve_status;
+    } cupdlpx_presolve_info_t;
+
+    cupdlpx_presolve_info_t *pslp_presolve(const lp_problem_t *original_prob, const pdhg_parameters_t *params);
+
+    cupdlpx_result_t *create_result_from_presolve(const cupdlpx_presolve_info_t *info, const lp_problem_t *original_prob);
+
+    const char *get_presolve_status_str(enum PresolveStatus_ status);
+
+    void 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,7 @@ 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 cupdlpx_result_t *result, const pdhg_parameters_t *params);
 
     void display_iteration_stats(const pdhg_solver_state_t *solver_state, bool verbose);
 
@@ -132,7 +129,7 @@ extern "C"
 
     void print_initial_feas_polish_info(bool is_primal_polish, const pdhg_parameters_t *params);
 
-    void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state, bool verbose,  bool is_primal_polish);
+    void display_feas_polish_iteration_stats(const pdhg_solver_state_t *state, bool verbose, bool is_primal_polish);
 
     void pdhg_feas_polish_final_log(const pdhg_solver_state_t *primal_state, const pdhg_solver_state_t *dual_state, bool verbose);
 

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "scikit_build_core.build"
 
 [project]
 name = "cupdlpx"
-version = "0.1.5"
+version = "0.2.0"
 description = "Python bindings for cuPDLPx (GPU-accelerated first-order LP solver)"
 readme = "README.md"
 license = { text = "Apache-2.0" }

python/README.md

@@ -153,10 +153,11 @@ Below is a list of commonly used parameters, their internal keys, and descriptio
 | `RestartNecessaryReduction` | `necessary_reduction_for_restart` | float | `0.5` | Necessary reduction factor required for a restart. |
 | `RestartKp` | `k_p` | float | `0.99` | Proportional coefficient for PID-controlled primal weight updates. |
 | `ReflectionCoeff` | `reflection_coefficient` | float | `1.0` | Reflection coefficient. |
-| `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.|
-| `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual.  |
 | `SVMaxIter` | `sv_max_iter` | int | 5000 | Maximum number of iterations for the power method |
 | `SVTol`| `sv_tol` | float | `1e-4` | Termination tolerance for the power method |
+| `Presolve`| `presolve` | float | `True` | Whether to use presolve. |
+| `FeasibilityPolishing` | `feasibility_polishing` | bool | `False` | Run feasibility polishing process.|
+| `FeasibilityPolishingTol` | `eps_feas_polish_relative` | float | `1e-6` | Relative tolerance for primal/dual residual.  |
 
 They can be set in multiple ways:
 

python/cupdlpx/PDLP.py

@@ -56,4 +56,6 @@
     # singular value estimation (power method)
     "SVMaxIter": "sv_max_iter",
     "SVTol": "sv_tol",
+    # presolve
+    "Presolve": "presolve",
 }

python_bindings/_core_bindings.cpp

@@ -170,6 +170,8 @@ static const char *status_to_str(termination_reason_t r)
         return "PRIMAL_INFEASIBLE";
     case TERMINATION_REASON_DUAL_INFEASIBLE:
         return "DUAL_INFEASIBLE";
+    case TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED:
+        return "INFEASIBLE_OR_UNBOUNDED";
     case TERMINATION_REASON_TIME_LIMIT:
         return "TIME_LIMIT";
     case TERMINATION_REASON_ITERATION_LIMIT:
@@ -198,6 +200,8 @@ static int status_to_code(termination_reason_t r)
         return 3;
     case TERMINATION_REASON_ITERATION_LIMIT:
         return 4;
+    case TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED:
+        return 5;
     case TERMINATION_REASON_UNSPECIFIED:
     default:
         return -1;
@@ -247,6 +251,9 @@ static py::dict get_default_params_py()
     d["sv_max_iter"] = p.sv_max_iter;
     d["sv_tol"] = p.sv_tol;
 
+    // presolve
+    d["presolve"] = p.presolve;
+
     return d;
 }
 
@@ -299,6 +306,9 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
     // power method for singular value estimation
     geti("sv_max_iter", p->sv_max_iter);
     getf("sv_tol", p->sv_tol);
+
+    // presolve
+    getb("presolve", p->presolve);
 }
 
 // view of matrix from Python

src/cli.c

@@ -16,6 +16,7 @@ limitations under the License.
 
 #include "cupdlpx.h"
 #include "mps_parser.h"
+#include "presolve.h"
 #include "solver.h"
 #include "utils.h"
 #include <getopt.h>
@@ -61,7 +62,7 @@ void save_solution(const double *data, int size, const char *output_dir,
                    const char *instance_name, const char *suffix)
 {
     char *file_path = get_output_path(output_dir, instance_name, suffix);
-    if (file_path == NULL)
+    if (file_path == NULL || data == NULL)
     {
         return;
     }
@@ -113,7 +114,19 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir,
     fprintf(outfile, "Absolute Objective Gap: %e\n", result->objective_gap);
     fprintf(outfile, "Relative Objective Gap: %e\n",
             result->relative_objective_gap);
-    if(result->feasibility_polishing_time > 0.0){
+    fprintf(outfile, "Rows: %d\n", result->num_constraints);
+    fprintf(outfile, "Columns: %d\n", result->num_variables);
+    fprintf(outfile, "Nonzeros: %d\n", result->num_nonzeros);
+    if (result->presolve_time > 0.0)
+    {
+        fprintf(outfile, "Presolve Status: %s\n", get_presolve_status_str(result->presolve_status));
+        fprintf(outfile, "Presolve Time (sec): %e\n", result->presolve_time);
+        fprintf(outfile, "Reduced Rows: %d\n", result->num_reduced_constraints);
+        fprintf(outfile, "Reduced Columns: %d\n", result->num_reduced_variables);
+        fprintf(outfile, "Reduced Nonzeros: %d\n", result->num_reduced_nonzeros);
+    }
+    if (result->feasibility_polishing_time > 0.0)
+    {
         fprintf(outfile, "Feasibility Polishing Time (sec): %e\n", result->feasibility_polishing_time);
         fprintf(outfile, "Feasibility Polishing Iteration Count: %d\n", result->feasibility_iteration);
     }
@@ -171,6 +184,8 @@ void print_usage(const char *prog_name)
                     "Enable feasibility use feasibility polishing (default: false).\n");
     fprintf(stderr, "      --eps_feas_polish <tolerance>   Relative feasibility "
                     "polish tolerance (default: 1e-6).\n");
+    fprintf(stderr, "      --no_presolve                   "
+                    "Disable presolve (default: enabled).\n");
 }
 
 int main(int argc, char *argv[])
@@ -195,10 +210,11 @@ int main(int argc, char *argv[])
         {"sv_max_iter", required_argument, 0, 1011},
         {"sv_tol", required_argument, 0, 1012},
         {"eval_freq", required_argument, 0, 1013},
+        {"no_presolve", no_argument, 0, 1014},
         {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)
         {
@@ -226,7 +242,7 @@ int main(int argc, char *argv[])
         case 1006: // --eps_feas_polish_relative
             params.termination_criteria.eps_feas_polish_relative = atof(optarg);
             break;
-        case 'f':  // --feasibility_polishing
+        case 'f': // --feasibility_polishing
             params.feasibility_polishing = true;
             break;
         case 1007: // --l_inf_ruiz_iter
@@ -250,6 +266,9 @@ int main(int argc, char *argv[])
         case 1013: // --eval_freq
             params.termination_evaluation_frequency = atoi(optarg);
             break;
+        case 1014: // --no_presolve
+            params.presolve = false;
+            break;
         case '?': // Unknown option
             return 1;
         }