diff --git a/.gitignore b/.gitignore index cd9fb58..866c33f 100644 --- a/.gitignore +++ b/.gitignore @@ -41,6 +41,7 @@ *.su *.idb *.pdb +Testing/* # Kernel Module Compile Results *.mod* @@ -60,3 +61,5 @@ test/* /.vscode /.venv /_b +*.whl +*.pyc diff --git a/README.md b/README.md index 320a9eb..daf779b 100644 --- a/README.md +++ b/README.md @@ -84,6 +84,7 @@ After building the project, the `./build/cupdlpx` binary can be invoked from the | `-v`, `--verbose` | `flag` | Enable verbose logging. | `false` | | `--time_limit` | `double` | Time limit in seconds. | `3600.0` | | `--iter_limit` | `int` | Iteration limit. | `2147483647` | +| `--opt_norm` | `string` | Norm for optimality criteria: `l2` or `linf` | `l2` | | `--eps_opt` | `double` | Relative optimality tolerance. | `1e-4` | | `--eps_feas` | `double` | Relative feasibility tolerance. | `1e-4` | | `--eps_infeas_detect` | `double` | Infeasibility detection tolerance. | `1e-10` | diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index b2cb101..c213b97 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -37,6 +37,12 @@ extern "C" TERMINATION_REASON_FEAS_POLISH_SUCCESS } termination_reason_t; + typedef enum + { + NORM_TYPE_L2 = 0, + NORM_TYPE_L_INF = 1 + } norm_type_t; + typedef struct { int num_variables; @@ -93,6 +99,7 @@ extern "C" restart_parameters_t restart_params; double reflection_coefficient; bool feasibility_polishing; + norm_type_t optimality_norm; bool presolve; } pdhg_parameters_t; diff --git a/internal/internal_types.h b/internal/internal_types.h index 1f42582..9eef900 100644 --- a/internal/internal_types.h +++ b/internal/internal_types.h @@ -103,6 +103,7 @@ typedef struct double initial_fixed_point_error; double last_trial_fixed_point_error; int inner_count; + norm_type_t optimality_norm; cusparseHandle_t sparse_handle; cublasHandle_t blas_handle; diff --git a/python/README.md b/python/README.md index 6e8e042..c5c1342 100644 --- a/python/README.md +++ b/python/README.md @@ -141,6 +141,7 @@ Below is a list of commonly used parameters, their internal keys, and descriptio | `IterationLimit` | `iteration_limit` | int | `2147483647` | Maximum number of iterations. | | `OutputFlag`, `LogToConsole` | `verbose` | bool | `False` | Enable (`True`) or disable (`False`) console logging output. | | `TermCheckFreq` | `termination_evaluation_frequency` | int | `200` | Frequency (in iterations) at which termination conditions are evaluated. | +| `OptimalityNorm` | `optimality_norm` | string | `"l2"` | Norm for optimality criteria. Use `"l2"` for L2 norm or `"linf"` for infinity norm. | | `OptimalityTol` | `eps_optimal_relative` | float | `1e-4` | Relative tolerance for optimality gap. Solver stops if the relative primal-dual gap ≤ this value. | | `FeasibilityTol` | `eps_feasible_relative` | float | `1e-4` | Relative feasibility tolerance for primal/dual residuals. | | `InfeasibleTol` | `eps_infeasible` | float | `1e-10` | Threshold for declaring infeasibility. | diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index 2a56350..0d2eec7 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -53,6 +53,8 @@ # feasibility polishing "FeasibilityPolishing": "feasibility_polishing", "FeasibilityPolishingTol": "eps_feas_polish_relative", + # termination criteria + "OptimalityNorm": "optimality_norm", # singular value estimation (power method) "SVMaxIter": "sv_max_iter", "SVTol": "sv_tol", diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index c22d66b..b18aba7 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -137,6 +137,21 @@ static const int32_t *get_index_ptr_i32(py::object obj, const char *name, throw std::invalid_argument(std::string(name) + " must be int32 or int64."); } +// helper function to convert norm string to enum +static norm_type_t parse_norm_string(const std::string& s) +{ + std::string lower = s; + std::transform(lower.begin(), lower.end(), lower.begin(), ::tolower); + + if (lower == "l2") { + return NORM_TYPE_L2; + } else if (lower == "linf") { + return NORM_TYPE_L_INF; + } else { + throw std::invalid_argument("Unknown norm type: " + s + ". Use 'l2' or 'linf'."); + } +} + // ensure 1D array or None with expected length static void ensure_len_or_null(py::object obj, const char *name, int expect_len) { @@ -247,6 +262,8 @@ static py::dict get_default_params_py() d["feasibility_polishing"] = p.feasibility_polishing; d["eps_feas_polish_relative"] = p.termination_criteria.eps_feas_polish_relative; + // Termination criteria norm + d["optimality_norm"] = (p.optimality_norm == NORM_TYPE_L_INF) ? "linf" : "l2"; // power method for singular value estimation d["sv_max_iter"] = p.sv_max_iter; d["sv_tol"] = p.sv_tol; @@ -270,6 +287,19 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p { if (d.contains(k)) tgt = py::cast(d[k]); }; auto getb = [&](const char *k, bool &tgt) { if (d.contains(k)) tgt = py::cast(d[k]); }; + auto get_norm = [&](const char *k, norm_type_t &tgt) + { + if (d.contains(k)) { + py::object val = d[k]; + if (py::isinstance(val)) { + std::string sval = py::cast(val); + tgt = parse_norm_string(sval); + } + else { + throw std::invalid_argument("optimality_norm must be a string ('l2'/'linf')"); + } + } + }; // verbosity getb("verbose", p->verbose); @@ -303,6 +333,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getb("feasibility_polishing", p->feasibility_polishing); getf("eps_feas_polish_relative", p->termination_criteria.eps_feas_polish_relative); + // Termination criteria norm + get_norm("optimality_norm", p->optimality_norm); // power method for singular value estimation geti("sv_max_iter", p->sv_max_iter); getf("sv_tol", p->sv_tol); diff --git a/src/cli.c b/src/cli.c index 690a9ff..f0bfe4e 100644 --- a/src/cli.c +++ b/src/cli.c @@ -201,6 +201,8 @@ void print_usage(const char *prog_name) "Enable feasibility use feasibility polishing (default: false).\n"); fprintf(stderr, " --eps_feas_polish Relative feasibility " "polish tolerance (default: 1e-6).\n"); + fprintf(stderr, " --opt_norm " + "Norm for optimality criteria: l2 or linf (default: l2).\n"); fprintf(stderr, " --no_presolve " "Disable presolve (default: enabled).\n"); } @@ -227,7 +229,8 @@ 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}, + {"opt_norm", required_argument, 0, 1014}, + {"no_presolve", no_argument, 0, 1015}, {0, 0, 0, 0}}; int opt; @@ -283,7 +286,20 @@ int main(int argc, char *argv[]) case 1013: // --eval_freq params.termination_evaluation_frequency = atoi(optarg); break; - case 1014: // --no_presolve + case 1014: // --opt_norm + { + const char *norm_str = optarg; + if (strcmp(norm_str, "l2") == 0) { + params.optimality_norm = NORM_TYPE_L2; + } else if (strcmp(norm_str, "linf") == 0) { + params.optimality_norm = NORM_TYPE_L_INF; + } else { + fprintf(stderr, "Error: opt_norm must be 'l2' or 'linf'\n"); + return 1; + } + } + break; + case 1015: // --no_presolve params.presolve = false; break; case '?': // Unknown option diff --git a/src/solver.cu b/src/solver.cu index 0e58dbf..f1abacd 100644 --- a/src/solver.cu +++ b/src/solver.cu @@ -224,6 +224,8 @@ initialize_solver_state(const pdhg_parameters_t *params, pdhg_solver_state_t *state = (pdhg_solver_state_t *)safe_calloc(1, sizeof(pdhg_solver_state_t)); + state->optimality_norm = params ? params->optimality_norm : NORM_TYPE_L2; + int n_vars = working_problem->num_variables; int n_cons = working_problem->num_constraints; size_t var_bytes = n_vars * sizeof(double); @@ -408,32 +410,59 @@ initialize_solver_state(const pdhg_parameters_t *params, free(temp_host); double sum_of_squares = 0.0; + double max_val = 0.0; + double val = 0.0; for (int i = 0; i < n_vars; ++i) { - sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i]; + if (state->optimality_norm == NORM_TYPE_L_INF) { + val = fabs(working_problem->objective_vector[i]); + if (val > max_val) max_val = val; + } else { + sum_of_squares += working_problem->objective_vector[i] * working_problem->objective_vector[i]; + } + } + + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->objective_vector_norm = max_val; + } else { + state->objective_vector_norm = sqrt(sum_of_squares); } - state->objective_vector_norm = sqrt(sum_of_squares); sum_of_squares = 0.0; + max_val = 0.0; + val = 0.0; for (int i = 0; i < n_cons; ++i) { double lower = working_problem->constraint_lower_bound[i]; double upper = working_problem->constraint_upper_bound[i]; - if (isfinite(lower) && (lower != upper)) - { - sum_of_squares += lower * lower; + if (state->optimality_norm == NORM_TYPE_L_INF) { + if (isfinite(lower) && (lower != upper)) { + val = fabs(lower); + if (val > max_val) max_val = val; + } + if (isfinite(upper)) { + val = fabs(upper); + if (val > max_val) max_val = val; + } + } else { + if (isfinite(lower) && (lower != upper)) { + sum_of_squares += lower * lower; + } + if (isfinite(upper)) { + sum_of_squares += upper * upper; + } } + } - if (isfinite(upper)) - { - sum_of_squares += upper * upper; - } + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->constraint_bound_norm = max_val; + } else { + state->constraint_bound_norm = sqrt(sum_of_squares); } - state->constraint_bound_norm = sqrt(sum_of_squares); state->num_blocks_primal = (state->num_variables + THREADS_PER_BLOCK - 1) / THREADS_PER_BLOCK; state->num_blocks_dual = diff --git a/src/utils.cu b/src/utils.cu index b86b307..0bd75d1 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -224,8 +224,8 @@ bool optimality_criteria_met(const pdhg_solver_state_t *state, double rel_opt_tol, double rel_feas_tol) { return state->relative_dual_residual < rel_feas_tol && - state->relative_primal_residual < rel_feas_tol && - state->relative_objective_gap < rel_opt_tol; + state->relative_primal_residual < rel_feas_tol && + state->relative_objective_gap < rel_opt_tol; } bool primal_infeasibility_criteria_met(const pdhg_solver_state_t *state, @@ -349,6 +349,7 @@ void set_default_parameters(pdhg_parameters_t *params) params->restart_params.k_d = 0.0; params->restart_params.i_smooth = 0.3; + params->optimality_norm = NORM_TYPE_L2; params->presolve = true; } @@ -707,13 +708,26 @@ void compute_residual(pdhg_solver_state_t *state) state->constraint_upper_bound_finite_val, state->num_constraints, state->num_variables); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, - state->primal_residual, 1, - &state->absolute_primal_residual)); + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, + state->num_constraints, state->primal_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, + state->primal_residual, 1, + &state->absolute_primal_residual)); + } + state->absolute_primal_residual /= state->constraint_bound_rescaling; - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, - state->dual_residual, 1, - &state->absolute_dual_residual)); + + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, + state->num_variables, state->dual_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, + state->dual_residual, 1, + &state->absolute_dual_residual)); + } + state->absolute_dual_residual /= state->objective_vector_rescaling; CUBLAS_CHECK(cublasDdot( @@ -736,12 +750,15 @@ void compute_residual(pdhg_solver_state_t *state) state->objective_vector_rescaling) + state->objective_constant; - state->relative_primal_residual = + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); + state->objective_gap = fabs(state->primal_objective_value - state->dual_objective_value); + state->relative_objective_gap = state->objective_gap / (1.0 + fabs(state->primal_objective_value) + fabs(state->dual_objective_value)); @@ -1230,8 +1247,17 @@ void compute_primal_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_ state->constraint_upper_bound, state->constraint_rescaling, state->num_constraints); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, state->primal_residual, 1, &state->absolute_primal_residual)); + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->absolute_primal_residual = get_vector_inf_norm(state->blas_handle, + state->num_constraints, state->primal_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_constraints, + state->primal_residual, 1, + &state->absolute_primal_residual)); + } + state->absolute_primal_residual /= state->constraint_bound_rescaling; + state->relative_primal_residual = state->absolute_primal_residual / (1.0 + state->constraint_bound_norm); CUBLAS_CHECK(cublasDdot(state->blas_handle, state->num_variables, ori_state->objective_vector, 1, state->pdhg_primal_solution, 1, &state->primal_objective_value)); @@ -1256,8 +1282,16 @@ void compute_dual_feas_polish_residual(pdhg_solver_state_t *state, const pdhg_so ori_state->constraint_upper_bound_finite_val, state->num_variables, state->num_constraints); - CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, state->dual_residual, 1, &state->absolute_dual_residual)); + if (state->optimality_norm == NORM_TYPE_L_INF) { + state->absolute_dual_residual = get_vector_inf_norm(state->blas_handle, + state->num_variables, state->dual_residual); + } else { + CUBLAS_CHECK(cublasDnrm2_v2_64(state->blas_handle, state->num_variables, + state->dual_residual, 1, &state->absolute_dual_residual)); + } + state->absolute_dual_residual /= state->objective_vector_rescaling; + state->relative_dual_residual = state->absolute_dual_residual / (1.0 + state->objective_vector_norm); double base_dual_objective; diff --git a/test/test.sh b/test/test.sh old mode 100644 new mode 100755