From 846822f9911a8bc78686a6a86c04ff511bbbbf14 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 24 Nov 2025 16:39:53 -0500 Subject: [PATCH 1/8] add PSLP presolver --- CMakeLists.txt | 12 +++ include/cupdlpx_types.h | 1 + internal/presolve.h | 31 ++++++++ src/cli.c | 9 ++- src/presolve.c | 167 ++++++++++++++++++++++++++++++++++++++++ src/solver.cu | 47 +++++++++-- 6 files changed, 261 insertions(+), 6 deletions(-) create mode 100644 internal/presolve.h create mode 100644 src/presolve.c diff --git a/CMakeLists.txt b/CMakeLists.txt index 517fbb5..75518a0 100644 --- a/CMakeLists.txt +++ b/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 ) # ----------------------------------------------------------------------------- diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 771c3e3..bb48653 100644 --- a/include/cupdlpx_types.h +++ b/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 diff --git a/internal/presolve.h b/internal/presolve.h new file mode 100644 index 0000000..c893901 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/cli.c b/src/cli.c index b9588ad..22546d6 100644 --- a/src/cli.c +++ b/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 #include #include @@ -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; } diff --git a/src/presolve.c b/src/presolve.c new file mode 100644 index 0000000..7f057bb --- /dev/null +++ b/src/presolve.c @@ -0,0 +1,167 @@ +#include "presolve.h" +#include "cupdlpx.h" +#include "PSLP_inf.h" +#include "PSLP_sol.h" +#include +#include +#include +#include +#include +#include + +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::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); +} \ No newline at end of file diff --git a/src/solver.cu b/src/solver.cu index a9e5f2b..f15ebdc 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -19,6 +19,7 @@ limitations under the License. #include "preconditioner.h" #include "solver.h" #include "utils.h" +#include "presolve.h" #include #include #include @@ -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 From ca4fda57a11deb7376566acea04783814f3a1f21 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Tue, 25 Nov 2025 17:52:07 -0500 Subject: [PATCH 2/8] fix: correct the reduced cost in PSLP postsolve --- include/cupdlpx_types.h | 1 + src/presolve.c | 12 +++--------- src/solver.cu | 42 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 46 insertions(+), 9 deletions(-) diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index bb48653..e50e1af 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -99,6 +99,7 @@ extern "C" double *primal_solution; double *dual_solution; + double *reduced_cost; int total_count; double rescaling_time_sec; diff --git a/src/presolve.c b/src/presolve.c index 7f057bb..440c7e8 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -44,8 +44,6 @@ lp_problem_t* convert_pslp_to_cupdlpx(PresolvedProblem *reduced_prob) { 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); @@ -58,7 +56,7 @@ lp_problem_t* convert_pslp_to_cupdlpx(PresolvedProblem *reduced_prob) { reduced_prob->rhs, // constraint upper bound reduced_prob->lbs, // variable lower bound reduced_prob->ubs, // variable upper bound - &obj_const + &reduced_prob->obj_offset ); return gpu_prob; @@ -81,6 +79,7 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const if (params->verbose) { info->settings->verbose = true; } + // info->settings->relax_bounds = false; // 3. Init Presolver info->presolver = new_presolver( @@ -113,7 +112,6 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const 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; @@ -133,16 +131,12 @@ cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 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->reduced_cost, 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; diff --git a/src/solver.cu b/src/solver.cu index f15ebdc..62eb798 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -204,6 +204,22 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, return final_result; } +__global__ void compute_and_rescale_reduced_cost_kernel( + double *reduced_cost, + const double *objective, + const double *dual_product, + const double *variable_rescaling, + const double objective_vector_rescaling, + const double constraint_bound_rescaling, + int n_vars) +{ + int i = blockIdx.x * blockDim.x + threadIdx.x; + if (i < n_vars) + { + reduced_cost[i] = (objective[i] - dual_product[i]) * variable_rescaling[i] / objective_vector_rescaling; + } +} + static pdhg_solver_state_t * initialize_solver_state(const lp_problem_t *original_problem, const rescale_info_t *rescale_info) @@ -951,12 +967,35 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state) cupdlpx_result_t *results = (cupdlpx_result_t *)safe_calloc(1, sizeof(cupdlpx_result_t)); + // Compute reduced cost + CUSPARSE_CHECK(cusparseDnVecSetValues(state->vec_dual_sol, + state->pdhg_dual_solution)); + CUSPARSE_CHECK(cusparseDnVecSetValues(state->vec_dual_prod, + state->dual_product)); + + CUSPARSE_CHECK(cusparseSpMV( + state->sparse_handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &HOST_ONE, + state->matAt, state->vec_dual_sol, &HOST_ZERO, state->vec_dual_prod, + CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG2, state->dual_spmv_buffer)); + + compute_and_rescale_reduced_cost_kernel<<num_blocks_primal, THREADS_PER_BLOCK>>>( + state->dual_slack, + state->objective_vector, + state->dual_product, + state->variable_rescaling, + state->objective_vector_rescaling, + state->constraint_bound_rescaling, + state->num_variables + ); + rescale_solution(state); results->primal_solution = (double *)safe_malloc(state->num_variables * sizeof(double)); results->dual_solution = (double *)safe_malloc(state->num_constraints * sizeof(double)); + results->reduced_cost = + (double *)safe_malloc(state->num_variables * sizeof(double)); CUDA_CHECK(cudaMemcpy(results->primal_solution, state->pdhg_primal_solution, state->num_variables * sizeof(double), @@ -964,6 +1003,9 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state) CUDA_CHECK(cudaMemcpy(results->dual_solution, state->pdhg_dual_solution, state->num_constraints * sizeof(double), cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(results->reduced_cost, state->dual_slack, + state->num_variables * sizeof(double), + cudaMemcpyDeviceToHost)); results->num_variables = state->num_variables; results->num_constraints = state->num_constraints; From 81e40bbbb87a1b1f4baaa3323ed74b54b254ffa7 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Tue, 25 Nov 2025 21:15:26 -0500 Subject: [PATCH 3/8] add presolve stats --- include/cupdlpx_types.h | 2 ++ internal/presolve.h | 3 +- internal/utils.h | 12 +++++--- src/cli.c | 23 +++++++++++++++ src/presolve.c | 13 +++++++-- src/solver.cu | 36 ++++++++++++++++++------ src/utils.cu | 62 +++++++++++++++++++++++++++++++++++------ 7 files changed, 127 insertions(+), 24 deletions(-) diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index e50e1af..b440601 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -16,6 +16,7 @@ limitations under the License. #pragma once +#include "PSLP_stats.h" #include #include @@ -120,6 +121,7 @@ extern "C" termination_reason_t termination_reason; double feasibility_polishing_time; int feasibility_iteration; + PresolveStats presolve_stats; } cupdlpx_result_t; // matrix formats diff --git a/internal/presolve.h b/internal/presolve.h index c893901..68949c1 100644 --- a/internal/presolve.h +++ b/internal/presolve.h @@ -2,7 +2,8 @@ #define PRESOLVE_H #include "cupdlpx.h" -#include "PSLP_API.h" +#include "PSLP_API.h" +#include "PSLP_stats.h" #ifdef __cplusplus extern "C" { diff --git a/internal/utils.h b/internal/utils.h index 4998801..dd3f709 100644 --- a/internal/utils.h +++ b/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); diff --git a/src/cli.c b/src/cli.c index 22546d6..7e9d499 100644 --- a/src/cli.c +++ b/src/cli.c @@ -140,6 +140,29 @@ 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_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); } diff --git a/src/presolve.c b/src/presolve.c index 440c7e8..8d59aab 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -124,14 +124,18 @@ cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 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 + + 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; return res; } if (!reduced_result || !info->presolver) return NULL; - postsolve(info->presolver, reduced_result->primal_solution, reduced_result->dual_solution, @@ -139,8 +143,13 @@ cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 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)); diff --git a/src/solver.cu b/src/solver.cu index 62eb798..5b1bf20 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -66,7 +66,7 @@ static void compute_next_pdhg_dual_solution(pdhg_solver_state_t *state); static void halpern_update(pdhg_solver_state_t *state, double reflection_coefficient); static void rescale_solution(pdhg_solver_state_t *state); -static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state); +static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, PresolveStats *presolve_stats); static void perform_restart(pdhg_solver_state_t *state, const pdhg_parameters_t *params); static void @@ -177,7 +177,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, state->total_count++; } - pdhg_final_log(state, params->verbose, state->termination_reason); + // pdhg_final_log(state, params->verbose, state->termination_reason); if (params->feasibility_polishing && state->termination_reason != TERMINATION_REASON_DUAL_INFEASIBLE && @@ -186,20 +186,33 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, feasibility_polish(params, state); } - cupdlpx_result_t *current_result = create_result_from_state(state); - + PresolveStats *stats_ptr = NULL; + if (presolve_info != NULL && presolve_info->presolver != NULL) { + stats_ptr = presolve_info->presolver->stats; + } + cupdlpx_result_t *current_result = create_result_from_state(state, stats_ptr); + 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; } + PresolveStats *stats_to_log = NULL; + if (presolve_info != NULL && presolve_info->presolver != NULL) { + stats_to_log = presolve_info->presolver->stats; + } + // TODO: use results for pdhg_final_log instead of state. + pdhg_final_log(state, stats_to_log, params->verbose, state->termination_reason); + + if (params->use_presolve && presolve_info) { + cupdlpx_result_free(current_result); + cupdlpx_presolve_info_free(presolve_info); + } + pdhg_solver_state_free(state); return final_result; } @@ -962,7 +975,7 @@ void rescale_info_free(rescale_info_t *info) free(info); } -static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state) +static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state, PresolveStats *presolve_stats) { cupdlpx_result_t *results = (cupdlpx_result_t *)safe_calloc(1, sizeof(cupdlpx_result_t)); @@ -1025,6 +1038,11 @@ static cupdlpx_result_t *create_result_from_state(pdhg_solver_state_t *state) results->termination_reason = state->termination_reason; results->feasibility_polishing_time = state->feasibility_polishing_time; results->feasibility_iteration = state->feasibility_iteration; + if (presolve_stats != NULL) { + results->presolve_stats = *presolve_stats; + } else { + memset(&(results->presolve_stats), 0, sizeof(PresolveStats)); + } return results; } diff --git a/src/utils.cu b/src/utils.cu index 172e731..8c39419 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -364,7 +364,27 @@ void print_initial_info(const pdhg_parameters_t *params, "------------------\n"); } -void pdhg_final_log(const pdhg_solver_state_t *state, bool verbose, +// void pdhg_final_log(const pdhg_solver_state_t *state, bool verbose, +// termination_reason_t reason) +// { +// if (verbose) +// { +// printf("-------------------------------------------------------------------" +// "--------------------\n"); +// } +// printf("Solution Summary\n"); +// printf(" Status : %s\n", termination_reason_to_string(reason)); +// printf(" Iterations : %d\n", state->total_count - 1); +// printf(" Solve time : %.3g sec\n", state->cumulative_time_sec); +// printf(" Primal obj : %.10g\n", state->primal_objective_value); +// printf(" Dual obj : %.10g\n", state->dual_objective_value); +// printf(" Primal infeas : %.3e\n", state->relative_primal_residual); +// printf(" Dual infeas : %.3e\n", state->relative_dual_residual); +// } + +void pdhg_final_log(const pdhg_solver_state_t *state, + const PresolveStats *stats, + bool verbose, termination_reason_t reason) { if (verbose) @@ -373,13 +393,39 @@ void pdhg_final_log(const pdhg_solver_state_t *state, bool verbose, "--------------------\n"); } printf("Solution Summary\n"); - printf(" Status : %s\n", termination_reason_to_string(reason)); - printf(" Iterations : %d\n", state->total_count - 1); - printf(" Solve time : %.3g sec\n", state->cumulative_time_sec); - printf(" Primal obj : %.10g\n", state->primal_objective_value); - printf(" Dual obj : %.10g\n", state->dual_objective_value); - printf(" Primal infeas : %.3e\n", state->relative_primal_residual); - printf(" Dual infeas : %.3e\n", state->relative_dual_residual); + printf(" Status : %s\n", termination_reason_to_string(reason)); + printf(" Iterations : %d\n", state->total_count - 1); + printf(" Solve time : %.3g sec\n", state->cumulative_time_sec); + printf(" Primal obj : %.10g\n", state->primal_objective_value); + printf(" Dual obj : %.10g\n", state->dual_objective_value); + printf(" Primal infeas : %.3e\n", state->relative_primal_residual); + printf(" Dual infeas : %.3e\n", state->relative_dual_residual); + + if (stats != NULL && stats->n_rows_original > 0) { + printf("\nPresolve Summary\n"); + printf(" [Dimensions]\n"); + printf(" Original : %d rows, %d cols, %d nnz\n", + stats->n_rows_original, stats->n_cols_original, stats->nnz_original); + printf(" Reduced : %d rows, %d cols, %d nnz\n", + stats->n_rows_reduced, stats->n_cols_reduced, stats->nnz_reduced); + + printf(" [Reduction Details (NNZ Removed)]\n"); + printf(" Trivial : %d\n", stats->nnz_removed_trivial); + printf(" Fast : %d\n", stats->nnz_removed_fast); + printf(" Primal Propagation : %d\n", stats->nnz_removed_primal_propagation); + printf(" Parallel Rows : %d\n", stats->nnz_removed_parallel_rows); + printf(" Parallel Cols : %d\n", stats->nnz_removed_parallel_cols); + + printf(" [Timing]\n"); + printf(" Total Presolve : %.3g sec\n", stats->presolve_total_time); + printf(" Init : %.3g sec\n", stats->ps_time_init); + printf(" Fast : %.3g sec\n", stats->ps_time_fast); + printf(" Medium : %.3g sec\n", stats->ps_time_medium); + printf(" Primal Propagation : %.3g sec\n", stats->ps_time_primal_propagation); + printf(" Parallel Rows : %.3g sec\n", stats->ps_time_parallel_rows); + printf(" Parallel Cols : %.3g sec\n", stats->ps_time_parallel_cols); + printf(" Postsolve : %.3g sec\n", stats->ps_time_post_solve); + } } void display_iteration_stats(const pdhg_solver_state_t *state, bool verbose) From 6ee68cf46cd804a021852bcba409a03c057294f0 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Wed, 26 Nov 2025 01:18:21 -0500 Subject: [PATCH 4/8] test PSLP relax_bounds false --- src/presolve.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/presolve.c b/src/presolve.c index 8d59aab..c7b2326 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -79,7 +79,7 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const if (params->verbose) { info->settings->verbose = true; } - // info->settings->relax_bounds = false; + info->settings->relax_bounds = false; // 3. Init Presolver info->presolver = new_presolver( From 4abff5465c04a6dd5d7e4c70a3f78fc39c3d61d2 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Sat, 29 Nov 2025 18:58:14 -0500 Subject: [PATCH 5/8] remove PSLP_INF --- src/presolve.c | 39 +++------------------------------------ 1 file changed, 3 insertions(+), 36 deletions(-) diff --git a/src/presolve.c b/src/presolve.c index c7b2326..8f87e80 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -1,6 +1,5 @@ #include "presolve.h" #include "cupdlpx.h" -#include "PSLP_inf.h" #include "PSLP_sol.h" #include #include @@ -9,27 +8,6 @@ #include #include -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::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)); @@ -44,11 +22,6 @@ lp_problem_t* convert_pslp_to_cupdlpx(PresolvedProblem *reduced_prob) { desc.data.csr.col_ind = reduced_prob->Ai; desc.data.csr.vals = reduced_prob->Ax; - 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, @@ -68,20 +41,14 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const 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 + // 1. Init Settings info->settings = default_settings(); if (params->verbose) { info->settings->verbose = true; } info->settings->relax_bounds = false; - // 3. Init Presolver + // 2. Init Presolver info->presolver = new_presolver( original_prob->constraint_matrix_values, original_prob->constraint_matrix_col_indices, @@ -98,7 +65,7 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const true ); - // 4. Run Presolve + // 3. Run Presolve PresolveStatus status = run_presolver(info->presolver); if (status & INFEASIBLE || status & UNBOUNDED) { From 49a6e8f44428bcfc211403dcbcf4c8b809a04261 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Sun, 30 Nov 2025 20:21:26 -0500 Subject: [PATCH 6/8] add presolve_setup_time and presolve_time record, improve reduce prob memory --- include/cupdlpx_types.h | 2 ++ internal/presolve.h | 1 + src/cli.c | 2 ++ src/presolve.c | 49 +++++++++++++++++++++-------------------- src/solver.cu | 7 ++++++ 5 files changed, 37 insertions(+), 24 deletions(-) diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index b440601..b21ffe0 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -105,6 +105,8 @@ extern "C" 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; diff --git a/internal/presolve.h b/internal/presolve.h index 68949c1..b12de54 100644 --- a/internal/presolve.h +++ b/internal/presolve.h @@ -15,6 +15,7 @@ typedef struct { 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); diff --git a/src/cli.c b/src/cli.c index 7e9d499..4c87c1e 100644 --- a/src/cli.c +++ b/src/cli.c @@ -141,6 +141,8 @@ void save_solver_summary(const cupdlpx_result_t *result, const char *output_dir, fprintf(outfile, "Feasibility Polishing Iteration Count: %d\n", result->feasibility_iteration); } if (result->presolve_stats.n_cols_original > 0) { + fprintf(outfile, "Presolve Time (sec): %e\n", result->presolve_time); + fprintf(outfile, "Presolve Setup Time (sec): %e\n", result->presolve_setup_time); 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); diff --git a/src/presolve.c b/src/presolve.c index 8f87e80..7c26181 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -1,6 +1,7 @@ #include "presolve.h" #include "cupdlpx.h" #include "PSLP_sol.h" +#include "utils.h" #include #include #include @@ -9,30 +10,29 @@ #include 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; - - 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 - &reduced_prob->obj_offset - ); - return gpu_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) { @@ -64,6 +64,7 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const info->settings, true ); + info->presolve_setup_time = (double)(clock() - start_time) / CLOCKS_PER_SEC; // 3. Run Presolve PresolveStatus status = run_presolver(info->presolver); @@ -130,7 +131,7 @@ cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 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->reduced_problem) lp_problem_free(info->reduced_problem); if (info->presolver) free_presolver(info->presolver); if (info->settings) free_settings(info->settings); free(info); diff --git a/src/solver.cu b/src/solver.cu index 5b1bf20..2215309 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -191,6 +191,10 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, stats_ptr = presolve_info->presolver->stats; } cupdlpx_result_t *current_result = create_result_from_state(state, stats_ptr); + if (presolve_info != NULL) { + current_result->presolve_time = presolve_info->presolve_time; + current_result->presolve_setup_time = presolve_info->presolve_setup_time; + } cupdlpx_result_t *final_result = NULL; @@ -207,6 +211,9 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, } // TODO: use results for pdhg_final_log instead of state. pdhg_final_log(state, stats_to_log, params->verbose, state->termination_reason); + printf("Presolve time (sec): %.3g\n", current_result->presolve_time); + printf("Presolve setup time (sec): %.3g\n", current_result->presolve_setup_time); + // printf(" Solve time : %.3g sec\n", state->cumulative_time_sec); if (params->use_presolve && presolve_info) { cupdlpx_result_free(current_result); From 86736275b4fa537956d86a628b5df4dcf95f2dc5 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Sun, 30 Nov 2025 20:22:41 -0500 Subject: [PATCH 7/8] set relax_bounds to True (default) --- src/presolve.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/presolve.c b/src/presolve.c index 7c26181..5fad158 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -46,7 +46,7 @@ cupdlpx_presolve_info_t* pslp_presolve(const lp_problem_t *original_prob, const if (params->verbose) { info->settings->verbose = true; } - info->settings->relax_bounds = false; + // info->settings->relax_bounds = false; // 2. Init Presolver info->presolver = new_presolver( From c8ba05d22a05999a03be809cfa749ce6670af43e Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Mon, 1 Dec 2025 00:26:23 -0500 Subject: [PATCH 8/8] fix Segmentation fault when presolve detected infeasible --- src/cli.c | 4 +++- src/presolve.c | 10 ++++++++++ src/solver.cu | 4 ++++ 3 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/cli.c b/src/cli.c index 4c87c1e..39aecca 100644 --- a/src/cli.c +++ b/src/cli.c @@ -140,9 +140,11 @@ 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_stats.n_cols_original > 0) { + 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); diff --git a/src/presolve.c b/src/presolve.c index 5fad158..1846107 100644 --- a/src/presolve.c +++ b/src/presolve.c @@ -99,6 +99,16 @@ cupdlpx_result_t* pslp_postsolve(cupdlpx_presolve_info_t *info, 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; } diff --git a/src/solver.cu b/src/solver.cu index 2215309..ba32ec4 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -114,6 +114,10 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params, if (presolve_info->problem_solved_during_presolve) { cupdlpx_result_t *early_result = pslp_postsolve(presolve_info, NULL, original_problem); + if (presolve_info != NULL) { + early_result->presolve_time = presolve_info->presolve_time; + early_result->presolve_setup_time = presolve_info->presolve_setup_time; + } cupdlpx_presolve_info_free(presolve_info); return early_result; }