Update PSLP to 0.0.4

CMakeLists.txt

@@ -8,7 +8,7 @@ project(cupdlpx LANGUAGES C CXX CUDA)
 
 set(CUPDLPX_VERSION_MAJOR 0)
 set(CUPDLPX_VERSION_MINOR 2)
-set(CUPDLPX_VERSION_PATCH 2)
+set(CUPDLPX_VERSION_PATCH 3)
 
 set(CUPDLPX_VERSION "${CUPDLPX_VERSION_MAJOR}.${CUPDLPX_VERSION_MINOR}.${CUPDLPX_VERSION_PATCH}")
 add_compile_definitions(CUPDLPX_VERSION="${CUPDLPX_VERSION}")
@@ -77,7 +77,7 @@ endif()
 
 include(FetchContent)
 
-set(PSLP_VERSION_TAG "v0.0.3")
+set(PSLP_VERSION_TAG "v0.0.4")
 
 FetchContent_Declare(
   pslp

internal/presolve.h

@@ -26,7 +26,7 @@ extern "C"
 
     const char *get_presolve_status_str(enum PresolveStatus_ status);
 
-    void pslp_postsolve(cupdlpx_presolve_info_t *info,
+    void pslp_postsolve(const cupdlpx_presolve_info_t *info,
                         cupdlpx_result_t *reduced_result,
                         const lp_problem_t *original_prob);
 

src/presolve.c

@@ -100,24 +100,21 @@ cupdlpx_presolve_info_t *pslp_presolve(const lp_problem_t *original_prob, const
     {
         printf("  %-15s : %s\n", "status", get_presolve_status_str(status));
         printf("  %-15s : %.3g sec\n", "presolve time", info->presolve_time);
+        printf("  %-15s : %d rows, %d columns, %d nonzeros\n",
+                "reduced problem",
+                info->presolver->reduced_prob->m,
+                info->presolver->reduced_prob->n,
+                info->presolver->reduced_prob->nnz);
     }
 
-    if (status & INFEASIBLE || status & UNBNDORINFEAS)
+    if (status & INFEASIBLE || status & UNBNDORINFEAS || info->presolver->reduced_prob->n == 0)
     {
         info->problem_solved_during_presolve = true;
         info->reduced_problem = NULL;
     }
     else
     {
         info->problem_solved_during_presolve = false;
-        if (params->verbose)
-        {
-            printf("  %-15s : %d rows, %d columns, %d nonzeros\n",
-                   "reduced problem",
-                   info->presolver->reduced_prob->m,
-                   info->presolver->reduced_prob->n,
-                   info->presolver->reduced_prob->nnz);
-        }
         info->reduced_problem = convert_pslp_to_cupdlpx(info->presolver->reduced_prob);
     }
     return info;
@@ -127,6 +124,14 @@ cupdlpx_result_t *create_result_from_presolve(const cupdlpx_presolve_info_t *inf
 {
 
     cupdlpx_result_t *result = (cupdlpx_result_t *)safe_calloc(1, sizeof(cupdlpx_result_t));
+    result->num_variables = original_prob->num_variables;
+    result->num_constraints = original_prob->num_constraints;
+    result->num_nonzeros = original_prob->constraint_matrix_num_nonzeros;
+    result->num_reduced_variables = info->presolver->reduced_prob->n;
+    result->num_reduced_constraints = info->presolver->reduced_prob->m;
+    result->num_reduced_nonzeros = info->presolver->reduced_prob->nnz;
+    result->presolve_status = info->presolve_status;
+    result->presolve_time = info->presolve_time;
 
     if (info->presolve_status == INFEASIBLE)
     {
@@ -136,20 +141,17 @@ cupdlpx_result_t *create_result_from_presolve(const cupdlpx_presolve_info_t *inf
     {
         result->termination_reason = TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED;
     }
+    else if (info->presolver->reduced_prob->n == 0)
+    {
+        result->termination_reason = TERMINATION_REASON_OPTIMAL;
+        pslp_postsolve(info, result, original_prob);
+        return result;
+    }
     else
     {
         result->termination_reason = TERMINATION_REASON_UNSPECIFIED;
     }
-    result->num_variables = original_prob->num_variables;
-    result->num_constraints = original_prob->num_constraints;
-    result->num_nonzeros = original_prob->constraint_matrix_num_nonzeros;
-    result->num_reduced_variables = info->presolver->reduced_prob->n;
-    result->num_reduced_constraints = info->presolver->reduced_prob->m;
-    result->num_reduced_nonzeros = info->presolver->reduced_prob->nnz;
-    result->presolve_status = info->presolve_status;
-    result->presolve_time = info->presolve_time;
     // result->presolve_stats = *(info->presolver->stats);
-    // TODO: Verify if setting solution pointers to NULL affects Python/Julia bindings.
     if (result->num_variables > 0)
     {
         result->primal_solution = (double *)safe_calloc(result->num_variables, sizeof(double));
@@ -162,15 +164,14 @@ cupdlpx_result_t *create_result_from_presolve(const cupdlpx_presolve_info_t *inf
     return result;
 }
 
-void pslp_postsolve(cupdlpx_presolve_info_t *info,
+void pslp_postsolve(const cupdlpx_presolve_info_t *info,
                     cupdlpx_result_t *result,
                     const lp_problem_t *original_prob)
 {
     postsolve(info->presolver,
               result->primal_solution,
               result->dual_solution,
-              result->reduced_cost,
-              result->primal_objective_value);
+              result->reduced_cost);
 
     result->num_reduced_variables = info->presolver->reduced_prob->n;
     result->num_reduced_constraints = info->presolver->reduced_prob->m;
@@ -184,11 +185,22 @@ void pslp_postsolve(cupdlpx_presolve_info_t *info,
     memcpy(result->primal_solution, info->presolver->sol->x, original_prob->num_variables * sizeof(double));
     memcpy(result->dual_solution, info->presolver->sol->y, original_prob->num_constraints * sizeof(double));
     memcpy(result->reduced_cost, info->presolver->sol->z, original_prob->num_variables * sizeof(double));
-    // result->primal_objective_value = info->presolver->sol->obj; // This is a bug in PSLP. We don't need to updated primal_objective_value since offset has been updated during presolve. Therefore, the original problem and reduced problem have the same objective value.
     result->presolve_time = info->presolve_time;
+    if (info->presolver->reduced_prob->n == 0)
+    {
+        double obj = 0.0;
+        for (int i = 0; i < original_prob->num_variables; i++)
+        {
+            obj += original_prob->objective_vector[i] * result->primal_solution[i];
+        }
+        obj += original_prob->objective_constant;
+        result->primal_objective_value = obj;
+        result->dual_objective_value = obj;
+    }
     // if (info->presolver->stats != NULL) {
     //     result->presolve_stats = *(info->presolver->stats);
     // }
+    return;
 }
 
 void cupdlpx_presolve_info_free(cupdlpx_presolve_info_t *info)

src/solver.cu

@@ -197,6 +197,7 @@ cupdlpx_result_t *optimize(const pdhg_parameters_t *params,
 
     pdhg_final_log(result, params);
     pdhg_solver_state_free(state);
+    CUDA_CHECK(cudaGetLastError());
     return result;
 }