Add matrix_zero_tol

docs/C_API.md

@@ -19,7 +19,7 @@ lp_problem_t *create_lp_problem(
 );
 
 cupdlpx_result_t* solve_lp_problem(
-    const lp_problem_t* prob,
+    lp_problem_t* prob,
     const pdhg_parameters_t* params    // NULL → use default parameters
 );
 ```
@@ -35,7 +35,7 @@ cupdlpx_result_t* solve_lp_problem(
 
 
 `solve_lp_problem` parameters:
-- `prob`: An LP problem built with `create_LP_problem`.
+- `prob`: An LP problem built with `create_LP_problem`. The solver may clean up the matrix (e.g., drop near-zero entries), so the struct must be mutable.
 - `params`: Solver parameters. If `NULL`, the solver will use default parameters.
 
 #### Example: Solving a Small LP
@@ -59,7 +59,6 @@ int main() {
     matrix_desc_t A_desc;
     A_desc.m = m; A_desc.n = n;
     A_desc.fmt = matrix_dense;
-    A_desc.zero_tolerance = 0.0;
     A_desc.data.dense.A = &A[0][0];
 
     // Objective coefficients

include/cupdlpx.h

@@ -38,7 +38,7 @@ extern "C"
 
     // solve the LP problem using PDHG
     cupdlpx_result_t *solve_lp_problem(
-        const lp_problem_t *prob,
+        lp_problem_t *prob,
         const pdhg_parameters_t *params);
 
     // parameter

include/cupdlpx_types.h

@@ -100,6 +100,7 @@ extern "C"
 		bool feasibility_polishing;
 		norm_type_t optimality_norm;
 		bool presolve;
+		double matrix_zero_tol;
 	} pdhg_parameters_t;
 
 	typedef struct
@@ -156,9 +157,6 @@ extern "C"
 		int n; // num_variables
 		matrix_format_t fmt;
 
-		// treat abs(x) < zero_tolerance as zero
-		double zero_tolerance;
-
 		union MatrixData
 		{
 			struct MatrixDense

internal/solver.h

@@ -25,7 +25,7 @@ extern "C"
 
     cupdlpx_result_t *optimize(
         const pdhg_parameters_t *params,
-        const lp_problem_t *original_problem);
+        lp_problem_t *original_problem);
 
 #ifdef __cplusplus
 }

internal/utils.h

@@ -97,7 +97,7 @@ extern "C"
         pdhg_solver_state_t *solver_state,
         const termination_criteria_t *criteria);
 
-    void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *problem);
+    void print_initial_info(const pdhg_parameters_t *params, lp_problem_t *problem);
 
     void pdhg_final_log(const cupdlpx_result_t *result, const pdhg_parameters_t *params);
 
@@ -117,10 +117,10 @@ extern "C"
                      int **row_ptr, int **col_ind, double **vals, int *nnz_out);
 
     int csc_to_csr(const matrix_desc_t *desc,
-                   int **row_ptr, int **col_ind, double **vals, int *nnz_out);
+                   int **row_ptr, int **col_ind, double **vals);
 
     int coo_to_csr(const matrix_desc_t *desc,
-                   int **row_ptr, int **col_ind, double **vals, int *nnz_out);
+                   int **row_ptr, int **col_ind, double **vals);
 
     void check_feas_polishing_termination_criteria(
         pdhg_solver_state_t *solver_state,

python/cupdlpx/PDLP.py

@@ -59,4 +59,5 @@
     "SVTol": "sv_tol",
     # presolve
     "Presolve": "presolve",
+    "MatrixZeroTol": "matrix_zero_tol",
 }

python/cupdlpx/model.py

@@ -368,7 +368,6 @@ def optimize(self):
             self.ub,
             self.constr_lb,
             self.constr_ub,
-            zero_tolerance=0.0,
             params=self._params,
             primal_start=self._primal_start,
             dual_start=self._dual_start

python_bindings/_core_bindings.cpp

@@ -270,6 +270,8 @@ static py::dict get_default_params_py()
     // presolve
     d["presolve"] = p.presolve;
 
+    d["matrix_zero_tol"] = p.matrix_zero_tol;
+
     return d;
 }
 
@@ -339,15 +341,16 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p
 
     // presolve
     getb("presolve", p->presolve);
+
+    getf("matrix_zero_tol", p->matrix_zero_tol);
 }
 
 // view of matrix from Python
-static PyMatrixView get_matrix_from_python(py::object A, double zero_tol)
+static PyMatrixView get_matrix_from_python(py::object A)
 {
     // initialize output
     PyMatrixView out;
     auto &desc = out.desc;
-    desc.zero_tolerance = zero_tol;
     // get shape
     if (!py::hasattr(A, "shape"))
     {
@@ -441,14 +444,13 @@ static py::dict solve_once(
     py::object variable_upper_bound,      // ub (optional → inf)
     py::object constraint_lower_bound,    // l  (optional → -inf)
     py::object constraint_upper_bound,    // u  (optional → inf)
-    double zero_tolerance = 0.0,          // zero filter tolerance
     py::object params = py::none(),       // PDHG parameters (optional → default)
     py::object primal_start = py::none(), // warm start primal solution (optional)
     py::object dual_start = py::none()    // warm start dual solution (optional)
 )
 {
     // parse matrix
-    PyMatrixView view = get_matrix_from_python(A, zero_tolerance);
+    PyMatrixView view = get_matrix_from_python(A);
     const int m = view.desc.m;
     const int n = view.desc.n;
     // get vector pointers
@@ -570,7 +572,6 @@ PYBIND11_MODULE(_cupdlpx_core, m)
           py::arg("variable_upper_bound") = py::none(),
           py::arg("constraint_lower_bound") = py::none(),
           py::arg("constraint_upper_bound") = py::none(),
-          py::arg("zero_tolerance") = 0.0,
           py::arg("params") = py::none(),
           py::arg("primal_start") = py::none(),
           py::arg("dual_start") = py::none());

src/cli.c

@@ -204,6 +204,8 @@ void print_usage(const char *prog_name)
                     "Norm for optimality criteria: l2 or linf (default: l2).\n");
     fprintf(stderr, "      --no_presolve                   "
                     "Disable presolve (default: enabled).\n");
+    fprintf(stderr, "      --matrix_zero_tol <tolerance>.  "
+                    "Zero tolerance in constraint matrix.\n");
 }
 
 int main(int argc, char *argv[])
@@ -229,6 +231,7 @@ int main(int argc, char *argv[])
         {"eval_freq", required_argument, 0, 1013},
         {"opt_norm", required_argument, 0, 1014},
         {"no_presolve", no_argument, 0, 1015},
+        {"matrix_zero_tol", required_argument, 0, 1016},
         {0, 0, 0, 0}};
 
     int opt;
@@ -297,6 +300,9 @@ int main(int argc, char *argv[])
         case 1015: // --no_presolve
             params.presolve = false;
             break;
+        case 1016: // --matrix_zero_tol
+            params.matrix_zero_tol = atof(optarg);
+            break;
         case '?': // Unknown option
             return 1;
         }
@@ -329,7 +335,7 @@ int main(int argc, char *argv[])
         return 1;
     }
 
-    cupdlpx_result_t *result = optimize(&params, problem);
+    cupdlpx_result_t *result = solve_lp_problem(problem, &params);
 
     if (result == NULL)
     {

src/cupdlpx.c

@@ -50,14 +50,13 @@ lp_problem_t *create_lp_problem(const double *objective_c,
     {
         int *row_ptr = NULL, *col_ind = NULL;
         double *vals = NULL;
-        int nnz = 0;
-        if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0)
+        if (csc_to_csr(A_desc, &row_ptr, &col_ind, &vals) != 0)
         {
             fprintf(stderr, "[interface] CSC->CSR failed.\n");
             free(prob);
             return NULL;
         }
-        prob->constraint_matrix_num_nonzeros = nnz;
+        prob->constraint_matrix_num_nonzeros = A_desc->data.csc.nnz;
         prob->constraint_matrix_row_pointers = row_ptr;
         prob->constraint_matrix_col_indices = col_ind;
         prob->constraint_matrix_values = vals;
@@ -68,14 +67,13 @@ lp_problem_t *create_lp_problem(const double *objective_c,
     {
         int *row_ptr = NULL, *col_ind = NULL;
         double *vals = NULL;
-        int nnz = 0;
-        if (coo_to_csr(A_desc, &row_ptr, &col_ind, &vals, &nnz) != 0)
+        if (coo_to_csr(A_desc, &row_ptr, &col_ind, &vals) != 0)
         {
             fprintf(stderr, "[interface] COO->CSR failed.\n");
             free(prob);
             return NULL;
         }
-        prob->constraint_matrix_num_nonzeros = nnz;
+        prob->constraint_matrix_num_nonzeros = A_desc->data.coo.nnz;
         prob->constraint_matrix_row_pointers = row_ptr;
         prob->constraint_matrix_col_indices = col_ind;
         prob->constraint_matrix_values = vals;
@@ -184,7 +182,7 @@ void set_start_values(lp_problem_t *prob, const double *primal,
     }
 }
 
-cupdlpx_result_t *solve_lp_problem(const lp_problem_t *prob,
+cupdlpx_result_t *solve_lp_problem(lp_problem_t *prob,
                                    const pdhg_parameters_t *params)
 {
     // argument checks

src/solver.cu

@@ -125,10 +125,10 @@ static void sync_inner_count_to_gpu(pdhg_solver_state_t *state);
 static void check_params_validity(const pdhg_parameters_t *params);
 
 cupdlpx_result_t *optimize(const pdhg_parameters_t *params,
-                           const lp_problem_t *original_problem)
+                           lp_problem_t *original_problem)
 {
-    print_initial_info(params, original_problem);
     check_params_validity(params);
+    print_initial_info(params, original_problem);
 
     cupdlpx_presolve_info_t *presolve_info = NULL;
     const lp_problem_t *working_problem = original_problem;