From 15955de2044c8a66f95d0f1fb25424bda9b4bc81 Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Thu, 20 Aug 2026 13:59:09 -0400 Subject: [PATCH 1/3] fix: MPS reader integer default bounds --- src/mps_parser.c | 151 ++++++++++++++++++++++++++++++++++++------ test/test_read_mps.py | 114 +++++++++++++++++++++++++++++++ 2 files changed, 246 insertions(+), 19 deletions(-) diff --git a/src/mps_parser.c b/src/mps_parser.c index 5ae437b..e98fb7a 100644 --- a/src/mps_parser.c +++ b/src/mps_parser.c @@ -26,6 +26,7 @@ limitations under the License. #include #define READER_BUFFER_SIZE (4 * 1024 * 1024) +#define MPS_INFINITE_BOUND 1e20 typedef struct NameNode { @@ -319,6 +320,9 @@ typedef struct double *var_upper_bounds; double *constraint_lower_bounds; double *constraint_upper_bounds; + unsigned char *col_binary_default; + unsigned char *col_has_lower; + int in_integer_block; size_t col_capacity; size_t constraint_capacity; @@ -365,12 +369,17 @@ static bool ensure_column_capacity(MpsParserState *state) state->objective_coeffs = (double *)safe_realloc(state->objective_coeffs, new_cap * sizeof(double)); state->var_lower_bounds = (double *)safe_realloc(state->var_lower_bounds, new_cap * sizeof(double)); state->var_upper_bounds = (double *)safe_realloc(state->var_upper_bounds, new_cap * sizeof(double)); + state->col_binary_default = + (unsigned char *)safe_realloc(state->col_binary_default, new_cap * sizeof(unsigned char)); + state->col_has_lower = (unsigned char *)safe_realloc(state->col_has_lower, new_cap * sizeof(unsigned char)); for (size_t i = state->col_capacity; i < new_cap; ++i) { state->objective_coeffs[i] = 0.0; state->var_lower_bounds[i] = 0.0; state->var_upper_bounds[i] = INFINITY; + state->col_binary_default[i] = 0; + state->col_has_lower[i] = 0; } state->col_capacity = new_cap; @@ -395,9 +404,46 @@ typedef enum SEC_RANGES, SEC_BOUNDS, SEC_OBJSENSE, + SEC_SOS, + SEC_UNSUPPORTED, SEC_ENDATA } MpsSection; +/* Sections that change the model beyond an LP: refuse the file rather than + * silently dropping them (HiGHS does the same for the ones it cannot parse). */ +static const char *const MPS_UNSUPPORTED_SECTIONS[] = {"QUADOBJ", + "QMATRIX", + "QSECTION", + "QCMATRIX", + "CSECTION", + "INDICATORS", + "DELAYEDROWS", + "MODELCUTS", + "USERCUTS", + "GENCONS", + "PWLOBJ", + "PWLNAM", + "PWLCON"}; + +/* Bounds at or beyond +/-1e20 are the conventional way of writing infinity in + * MPS files; HiGHS applies the same threshold (options.infinite_bound). */ +static int apply_infinite_bounds(double *lower, double *upper, size_t n, const char *what) +{ + for (size_t i = 0; i < n; ++i) + { + if (lower[i] >= MPS_INFINITE_BOUND || upper[i] <= -MPS_INFINITE_BOUND) + { + fprintf(stderr, "ERROR: %s %zu has bounds [%g, %g].\n", what, i, lower[i], upper[i]); + return -1; + } + if (lower[i] <= -MPS_INFINITE_BOUND) + lower[i] = -INFINITY; + if (upper[i] >= MPS_INFINITE_BOUND) + upper[i] = INFINITY; + } + return 0; +} + lp_problem_t *read_mps_file(const char *filename) { MpsParserState state = {0}; @@ -437,22 +483,42 @@ lp_problem_t *read_mps_file(const char *filename) if (isalpha((unsigned char)tokens[0][0])) { + /* A section header is alone on its line; OBJSENSE may carry its value. + * Checking the token count first keeps this off the hot path. */ MpsSection next_section = SEC_NONE; - if (strcmp(tokens[0], "ROWS") == 0) - next_section = SEC_ROWS; - else if (strcmp(tokens[0], "COLUMNS") == 0) - next_section = SEC_COLUMNS; - else if (strcmp(tokens[0], "RHS") == 0) - next_section = SEC_RHS; - else if (strcmp(tokens[0], "RANGES") == 0) - next_section = SEC_RANGES; - else if (strcmp(tokens[0], "BOUNDS") == 0) - next_section = SEC_BOUNDS; + if (n_tokens == 1) + { + if (strcmp(tokens[0], "ROWS") == 0) + next_section = SEC_ROWS; + else if (strcmp(tokens[0], "COLUMNS") == 0) + next_section = SEC_COLUMNS; + else if (strcmp(tokens[0], "RHS") == 0) + next_section = SEC_RHS; + else if (strcmp(tokens[0], "RANGES") == 0) + next_section = SEC_RANGES; + else if (strcmp(tokens[0], "BOUNDS") == 0) + next_section = SEC_BOUNDS; + else if (strcmp(tokens[0], "OBJSENSE") == 0 || strcmp(tokens[0], "OBJSENS") == 0) + next_section = SEC_OBJSENSE; + else if (strcmp(tokens[0], "SOS") == 0 || strcmp(tokens[0], "SETS") == 0) + next_section = SEC_SOS; + else if (strcmp(tokens[0], "ENDATA") == 0) + next_section = SEC_ENDATA; + else + { + for (size_t k = 0; k < sizeof(MPS_UNSUPPORTED_SECTIONS) / sizeof(MPS_UNSUPPORTED_SECTIONS[0]); ++k) + { + if (strcmp(tokens[0], MPS_UNSUPPORTED_SECTIONS[k]) == 0) + { + next_section = SEC_UNSUPPORTED; + break; + } + } + } + } else if (strcmp(tokens[0], "OBJSENSE") == 0 || strcmp(tokens[0], "OBJSENS") == 0) - next_section = SEC_OBJSENSE; - else if (strcmp(tokens[0], "ENDATA") == 0) { - next_section = SEC_ENDATA; + next_section = SEC_OBJSENSE; } bool inline_max = next_section == SEC_OBJSENSE && n_tokens >= 2 && @@ -463,6 +529,12 @@ lp_problem_t *read_mps_file(const char *filename) if (is_header) { + if (next_section == SEC_UNSUPPORTED) + { + fprintf(stderr, "ERROR: MPS file reader cannot parse %s section.\n", tokens[0]); + state.error_flag = 1; + break; + } if (current_section == SEC_ROWS && next_section != SEC_ROWS && !rows_finalized) { if (finalize_rows(&state) != 0) @@ -515,6 +587,8 @@ lp_problem_t *read_mps_file(const char *filename) if (parse_bounds_section(&state, tokens, n_tokens) != 0) state.error_flag = 1; break; + case SEC_SOS: + break; default: break; @@ -530,6 +604,23 @@ lp_problem_t *read_mps_file(const char *filename) return NULL; } + for (size_t i = 0; i < state.col_map.size; ++i) + { + if (state.col_binary_default[i]) + { + state.var_lower_bounds[i] = 0.0; + state.var_upper_bounds[i] = 1.0; + } + } + + if (apply_infinite_bounds(state.var_lower_bounds, state.var_upper_bounds, state.col_map.size, "Column") != 0 || + apply_infinite_bounds( + state.constraint_lower_bounds, state.constraint_upper_bounds, state.row_map.size, "Row") != 0) + { + free_parser_state(&state); + return NULL; + } + lp_problem_t *prob = safe_calloc(1, sizeof(lp_problem_t)); prob->num_variables = state.col_map.size; @@ -600,10 +691,8 @@ static int finalize_rows(MpsParserState *state) } } - if (obj_idx == -1 && state->num_buffered_rows > 0) - { - obj_idx = 0; - } + if (obj_idx == -1) + fprintf(stderr, "WARNING: No objective (N) row found in ROWS section; objective is zero.\n"); if (obj_idx != -1) { @@ -667,6 +756,13 @@ static int parse_columns_section(MpsParserState *state, char **tokens, int n_tok if (n_tokens >= 2 && strcmp(tokens[1], "'MARKER'") == 0) { + if (n_tokens >= 3) + { + if (strcmp(tokens[2], "'INTORG'") == 0) + state->in_integer_block = 1; + else if (strcmp(tokens[2], "'INTEND'") == 0) + state->in_integer_block = 0; + } return 0; } @@ -697,9 +793,12 @@ static int parse_columns_section(MpsParserState *state, char **tokens, int n_tok if (!ensure_column_capacity(state)) return -1; + size_t n_cols_before = state->col_map.size; int col_idx = namemap_put(&state->col_map, col_name); if (col_idx == -1) return -1; + if (state->col_map.size > n_cols_before && state->in_integer_block) + state->col_binary_default[col_idx] = 1; for (int i = pair_start_index; i + 1 < n_tokens; i += 2) { @@ -811,27 +910,38 @@ static int parse_bounds_section(MpsParserState *state, char **tokens, int n_toke if (col_idx == -1) return 0; - if (strcmp(bound_type, "LO") == 0) + /* Any BOUNDS entry cancels the implicit [0, 1] default of an integer column. */ + state->col_binary_default[col_idx] = 0; + + if (strcmp(bound_type, "LO") == 0 || strcmp(bound_type, "LI") == 0) { state->var_lower_bounds[col_idx] = value; + state->col_has_lower[col_idx] = 1; } - else if (strcmp(bound_type, "UP") == 0) + else if (strcmp(bound_type, "UP") == 0 || strcmp(bound_type, "UI") == 0) { state->var_upper_bounds[col_idx] = value; + /* A negative upper bound on a column with no explicit lower bound means + * the lower bound is -infinity (CPLEX/Gurobi/SCIP convention). */ + if (value < 0.0 && !state->col_has_lower[col_idx]) + state->var_lower_bounds[col_idx] = -INFINITY; } else if (strcmp(bound_type, "FX") == 0) { state->var_lower_bounds[col_idx] = value; state->var_upper_bounds[col_idx] = value; + state->col_has_lower[col_idx] = 1; } else if (strcmp(bound_type, "FR") == 0) { state->var_lower_bounds[col_idx] = -INFINITY; state->var_upper_bounds[col_idx] = INFINITY; + state->col_has_lower[col_idx] = 1; } else if (strcmp(bound_type, "MI") == 0) { state->var_lower_bounds[col_idx] = -INFINITY; + state->col_has_lower[col_idx] = 1; } else if (strcmp(bound_type, "PL") == 0) { @@ -841,6 +951,7 @@ static int parse_bounds_section(MpsParserState *state, char **tokens, int n_toke { state->var_lower_bounds[col_idx] = 0.0; state->var_upper_bounds[col_idx] = 1.0; + state->col_has_lower[col_idx] = 1; } return 0; } @@ -905,6 +1016,8 @@ static void free_parser_state(MpsParserState *state) free(state->objective_coeffs); free(state->var_lower_bounds); free(state->var_upper_bounds); + free(state->col_binary_default); + free(state->col_has_lower); free(state->constraint_lower_bounds); free(state->constraint_upper_bounds); free(state->objective_row_name); diff --git a/test/test_read_mps.py b/test/test_read_mps.py index 6167fa6..827c6f5 100644 --- a/test/test_read_mps.py +++ b/test/test_read_mps.py @@ -162,3 +162,117 @@ def test_read_and_optimize_maximize(mps_max_file, atol): assert model.Status == PDLP.OPTIMAL, f"Unexpected termination status: {model.Status}" assert np.allclose(model.X, [1.5, 1.75], atol=atol), f"Unexpected primal solution: {model.X}" assert np.isclose(model.ObjVal, 3.25, atol=atol), f"Unexpected objective value: {model.ObjVal}" + + +# --------------------------------------------------------------------------- +# BOUNDS-section semantics +# --------------------------------------------------------------------------- + +MPS_BOUNDS = """NAME BOUNDS +ROWS + N COST + L R1 +COLUMNS + C0 COST 1.0 R1 1.0 + MARKER 'MARKER' 'INTORG' + I1 COST 1.0 R1 1.0 + I2 COST 1.0 R1 1.0 + I3 COST 1.0 R1 1.0 + I4 COST 1.0 R1 1.0 + I5 COST 1.0 R1 1.0 + MARKER 'MARKER' 'INTEND' + C6 COST 1.0 R1 1.0 + C7 COST 1.0 R1 1.0 + C8 COST 1.0 R1 1.0 +RHS + RHS R1 5.0 +BOUNDS + UP BND I2 5.0 + LO BND I3 2.0 + MI BND I4 + UP BND I5 -3.0 + UP BND C6 -3.0 + LO BND C7 -10.0 + UP BND C7 -3.0 + LO BND C8 -1e30 + UP BND C8 1e20 +ENDATA +""" + + +def _read_text(tmp_path, name, text): + path = tmp_path / name + path.write_text(text) + return cupdlpx.read(path) + + +def test_read_bounds_conventions(tmp_path): + """ + - integer (MARKER) columns with no BOUNDS entry default to [0, 1] + - any BOUNDS entry on such a column cancels that default + - negative UP with no explicit lower bound implies lb = -inf, but an + explicit lower bound (C7) is kept + - |bound| >= 1e20 is treated as infinite + """ + model = _read_text(tmp_path, "bounds.mps", MPS_BOUNDS) + inf = np.inf + # C0 I1 I2 I3 I4 I5 C6 C7 C8 + expected_lb = [0.0, 0.0, 0.0, 2.0, -inf, -inf, -inf, -10.0, -inf] + expected_ub = [inf, 1.0, 5.0, inf, inf, -3.0, -3.0, -3.0, inf] + assert np.array_equal(model.lb, expected_lb), model.lb + assert np.array_equal(model.ub, expected_ub), model.ub + + +MPS_NO_OBJ_ROW = """NAME NOOBJ +ROWS + E R1 + L R2 +COLUMNS + X1 R1 1.0 R2 1.0 + X2 R1 2.0 +RHS + RHS R1 5.0 R2 2.0 +ENDATA +""" + + +def test_read_without_objective_row(tmp_path): + """ + A file with no N row has a zero objective; all rows stay constraints + (previously the first row was silently taken as the objective). + """ + model = _read_text(tmp_path, "noobj.mps", MPS_NO_OBJ_ROW) + assert model.num_vars == 2 + assert model.num_constrs == 2 + assert np.allclose(model.c, [0.0, 0.0]) + assert np.allclose(model.constr_lb, [5.0, -np.inf]) + assert np.allclose(model.constr_ub, [5.0, 2.0]) + + +MPS_SOS = MPS_MIN.replace( + "ENDATA\n", + "SOS\n S1 SOS s1\n X1 1\n X2 2\nENDATA\n", +) + +MPS_QUADOBJ = MPS_MIN.replace( + "ENDATA\n", + "QUADOBJ\n X1 X1 2.0\nENDATA\n", +) + + +def test_read_sos_section_is_ignored(tmp_path): + """ + SOS only restricts the integer feasible set, so the LP relaxation is + unchanged and the section is skipped. + """ + model = _read_text(tmp_path, "sos.mps", MPS_SOS) + _check_min_model_data(model) + + +def test_read_unsupported_section_raises(tmp_path): + """ + Sections that change the model beyond an LP (quadratic terms, cones, + indicators, ...) must not be silently dropped. + """ + with pytest.raises(Exception): + _read_text(tmp_path, "quadobj.mps", MPS_QUADOBJ) From 69ebf74c05c6fa745cec2e93b43b42171714b1bf Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 21 Aug 2026 20:05:02 -0400 Subject: [PATCH 2/3] fix: MPS reader integer default bounds --- src/mps_parser.c | 28 ---------------------------- test/test_read_mps.py | 8 +++++--- 2 files changed, 5 insertions(+), 31 deletions(-) diff --git a/src/mps_parser.c b/src/mps_parser.c index e98fb7a..6119429 100644 --- a/src/mps_parser.c +++ b/src/mps_parser.c @@ -26,7 +26,6 @@ limitations under the License. #include #define READER_BUFFER_SIZE (4 * 1024 * 1024) -#define MPS_INFINITE_BOUND 1e20 typedef struct NameNode { @@ -425,25 +424,6 @@ static const char *const MPS_UNSUPPORTED_SECTIONS[] = {"QUADOBJ", "PWLNAM", "PWLCON"}; -/* Bounds at or beyond +/-1e20 are the conventional way of writing infinity in - * MPS files; HiGHS applies the same threshold (options.infinite_bound). */ -static int apply_infinite_bounds(double *lower, double *upper, size_t n, const char *what) -{ - for (size_t i = 0; i < n; ++i) - { - if (lower[i] >= MPS_INFINITE_BOUND || upper[i] <= -MPS_INFINITE_BOUND) - { - fprintf(stderr, "ERROR: %s %zu has bounds [%g, %g].\n", what, i, lower[i], upper[i]); - return -1; - } - if (lower[i] <= -MPS_INFINITE_BOUND) - lower[i] = -INFINITY; - if (upper[i] >= MPS_INFINITE_BOUND) - upper[i] = INFINITY; - } - return 0; -} - lp_problem_t *read_mps_file(const char *filename) { MpsParserState state = {0}; @@ -613,14 +593,6 @@ lp_problem_t *read_mps_file(const char *filename) } } - if (apply_infinite_bounds(state.var_lower_bounds, state.var_upper_bounds, state.col_map.size, "Column") != 0 || - apply_infinite_bounds( - state.constraint_lower_bounds, state.constraint_upper_bounds, state.row_map.size, "Row") != 0) - { - free_parser_state(&state); - return NULL; - } - lp_problem_t *prob = safe_calloc(1, sizeof(lp_problem_t)); prob->num_variables = state.col_map.size; diff --git a/test/test_read_mps.py b/test/test_read_mps.py index 827c6f5..aa003cb 100644 --- a/test/test_read_mps.py +++ b/test/test_read_mps.py @@ -212,13 +212,15 @@ def test_read_bounds_conventions(tmp_path): - any BOUNDS entry on such a column cancels that default - negative UP with no explicit lower bound implies lb = -inf, but an explicit lower bound (C7) is kept - - |bound| >= 1e20 is treated as infinite + + Large finite bounds (C8) are read as written; turning them into infinities + is the solver's infinite_bound parameter, not the reader's job. """ model = _read_text(tmp_path, "bounds.mps", MPS_BOUNDS) inf = np.inf # C0 I1 I2 I3 I4 I5 C6 C7 C8 - expected_lb = [0.0, 0.0, 0.0, 2.0, -inf, -inf, -inf, -10.0, -inf] - expected_ub = [inf, 1.0, 5.0, inf, inf, -3.0, -3.0, -3.0, inf] + expected_lb = [0.0, 0.0, 0.0, 2.0, -inf, -inf, -inf, -10.0, -1e30] + expected_ub = [inf, 1.0, 5.0, inf, inf, -3.0, -3.0, -3.0, 1e20] assert np.array_equal(model.lb, expected_lb), model.lb assert np.array_equal(model.ub, expected_ub), model.ub From 27ca1f2b46b2c419c9ee3a277939f5a0ec3ad9bd Mon Sep 17 00:00:00 2001 From: Zedong Peng Date: Fri, 21 Aug 2026 20:07:23 -0400 Subject: [PATCH 3/3] fix: treat large bounds as infinite and warn on huge coefficients --- include/cupdlpx_types.h | 1 + python/cupdlpx/PDLP.py | 1 + python/cupdlpx/model.py | 2 + python_bindings/_core_bindings.cpp | 4 ++ src/cli.c | 7 +++ src/utils.cu | 92 +++++++++++++++++++++++++++++- 6 files changed, 106 insertions(+), 1 deletion(-) diff --git a/include/cupdlpx_types.h b/include/cupdlpx_types.h index 2d6e148..9a014e9 100644 --- a/include/cupdlpx_types.h +++ b/include/cupdlpx_types.h @@ -108,6 +108,7 @@ extern "C" norm_type_t optimality_norm; bool presolve; double matrix_zero_tol; + double infinite_bound; } pdhg_parameters_t; typedef struct diff --git a/python/cupdlpx/PDLP.py b/python/cupdlpx/PDLP.py index e990adb..72b1e44 100644 --- a/python/cupdlpx/PDLP.py +++ b/python/cupdlpx/PDLP.py @@ -64,4 +64,5 @@ # presolve "Presolve": "presolve", "MatrixZeroTol": "matrix_zero_tol", + "InfiniteBound": "infinite_bound", } diff --git a/python/cupdlpx/model.py b/python/cupdlpx/model.py index 82d2a7e..ed8d700 100644 --- a/python/cupdlpx/model.py +++ b/python/cupdlpx/model.py @@ -62,6 +62,7 @@ "eps_feas_polish_relative", "sv_tol", "matrix_zero_tol", + "infinite_bound", } ) _POSITIVE_FLOAT_PARAMS = frozenset( @@ -70,6 +71,7 @@ "eps_feasible_relative", "eps_feas_polish_relative", "sv_tol", + "infinite_bound", } ) _NONNEGATIVE_FLOAT_PARAMS = frozenset({"time_sec_limit", "matrix_zero_tol"}) diff --git a/python_bindings/_core_bindings.cpp b/python_bindings/_core_bindings.cpp index 40971e7..16aa450 100644 --- a/python_bindings/_core_bindings.cpp +++ b/python_bindings/_core_bindings.cpp @@ -307,6 +307,7 @@ static py::dict get_default_params_py() d["presolve"] = p.presolve; d["matrix_zero_tol"] = p.matrix_zero_tol; + d["infinite_bound"] = p.infinite_bound; return d; } @@ -416,6 +417,7 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p getb("presolve", p->presolve); getf("matrix_zero_tol", p->matrix_zero_tol); + getf("infinite_bound", p->infinite_bound); if (p->termination_evaluation_frequency <= 0) throw std::invalid_argument("termination_evaluation_frequency must be positive."); @@ -435,6 +437,8 @@ static void parse_params_from_python(py::object params_obj, pdhg_parameters_t *p throw std::invalid_argument("sv_tol must be positive."); if (p->termination_criteria.time_sec_limit < 0.0) throw std::invalid_argument("time_sec_limit must be nonnegative."); + if (p->infinite_bound <= 0.0) + throw std::invalid_argument("infinite_bound must be positive."); if (p->matrix_zero_tol < 0.0) throw std::invalid_argument("matrix_zero_tol must be nonnegative."); } diff --git a/src/cli.c b/src/cli.c index d8491c0..4cb4f3f 100644 --- a/src/cli.c +++ b/src/cli.c @@ -212,6 +212,9 @@ void print_usage(const char *prog_name) fprintf(stderr, " --matrix_zero_tol . " "Zero tolerance in constraint matrix.\n"); + fprintf(stderr, + " --infinite_bound . " + "Bounds at or beyond this are treated as infinite (default: 1e20).\n"); } int main(int argc, char *argv[]) @@ -238,6 +241,7 @@ int main(int argc, char *argv[]) {"opt_norm", required_argument, 0, 1014}, {"no_presolve", no_argument, 0, 1015}, {"matrix_zero_tol", required_argument, 0, 1016}, + {"infinite_bound", required_argument, 0, 1017}, {0, 0, 0, 0}}; int opt; @@ -317,6 +321,9 @@ int main(int argc, char *argv[]) case 1016: // --matrix_zero_tol params.matrix_zero_tol = atof(optarg); break; + case 1017: // --infinite_bound + params.infinite_bound = atof(optarg); + break; case '?': // Unknown option return 1; } diff --git a/src/utils.cu b/src/utils.cu index 2f450e4..1e0d47c 100644 --- a/src/utils.cu +++ b/src/utils.cu @@ -329,8 +329,12 @@ void set_default_parameters(pdhg_parameters_t *params) params->optimality_norm = NORM_TYPE_L2; params->presolve = true; params->matrix_zero_tol = 1e-9; + params->infinite_bound = 1e20; } +#define MATRIX_LARGE_VALUE 1e15 +#define OBJECTIVE_LARGE_VALUE 1e20 + void filter_constraint_matrix_entries(lp_problem_t *out, const lp_problem_t *in, const pdhg_parameters_t *params) { if (out == NULL || in == NULL) @@ -358,17 +362,36 @@ void filter_constraint_matrix_entries(lp_problem_t *out, const lp_problem_t *in, } int filtered_nnz = 0; + int num_large = 0; + double max_large = 0.0; for (int i = 0; i < num_rows; ++i) { for (int k = row_ptr[i]; k < row_ptr[i + 1]; ++k) { - if (fabs(vals[k]) > params->matrix_zero_tol) + const double abs_value = fabs(vals[k]); + if (abs_value > params->matrix_zero_tol) { ++filtered_nnz; } + if (abs_value >= MATRIX_LARGE_VALUE) + { + ++num_large; + if (abs_value > max_large) + max_large = abs_value; + } } } + if (num_large > 0) + { + fprintf(stderr, + "WARNING: %d constraint matrix %s |value| >= %.1e (largest %.3e); the problem is badly scaled.\n", + num_large, + (num_large == 1 ? "entry has" : "entries have"), + MATRIX_LARGE_VALUE, + max_large); + } + if (filtered_nnz == nnz) { return; @@ -412,6 +435,27 @@ void filter_constraint_matrix_entries(lp_problem_t *out, const lp_problem_t *in, out->constraint_matrix_num_nonzeros = filtered_nnz; } +/* Bounds at or beyond +/-infinite_bound stand for an infinite bound. */ +static void replace_large_bounds_with_infinity( + const double *in_lower, const double *in_upper, double **lower, double **upper, int n, double infinite_bound) +{ + bool any = false; + for (int i = 0; i < n && !any; ++i) + any = in_lower[i] <= -infinite_bound || in_upper[i] >= infinite_bound; + if (!any) + return; + + double *new_lower = (double *)safe_malloc((size_t)n * sizeof(double)); + double *new_upper = (double *)safe_malloc((size_t)n * sizeof(double)); + for (int i = 0; i < n; ++i) + { + new_lower[i] = (in_lower[i] <= -infinite_bound) ? -INFINITY : in_lower[i]; + new_upper[i] = (in_upper[i] >= infinite_bound) ? INFINITY : in_upper[i]; + } + *lower = new_lower; + *upper = new_upper; +} + lp_problem_t preprocess_problem(const lp_problem_t *original, const pdhg_parameters_t *params) { lp_problem_t working = *original; @@ -426,6 +470,41 @@ lp_problem_t preprocess_problem(const lp_problem_t *original, const pdhg_paramet working.objective_constant = -original->objective_constant; working.objective_sense = OBJECTIVE_SENSE_MINIMIZE; } + replace_large_bounds_with_infinity(original->variable_lower_bound, + original->variable_upper_bound, + &working.variable_lower_bound, + &working.variable_upper_bound, + original->num_variables, + params->infinite_bound); + replace_large_bounds_with_infinity(original->constraint_lower_bound, + original->constraint_upper_bound, + &working.constraint_lower_bound, + &working.constraint_upper_bound, + original->num_constraints, + params->infinite_bound); + + int num_large_obj = 0; + double max_large_obj = 0.0; + for (int i = 0; i < original->num_variables; ++i) + { + const double abs_cost = fabs(original->objective_vector[i]); + if (abs_cost >= OBJECTIVE_LARGE_VALUE) + { + ++num_large_obj; + if (abs_cost > max_large_obj) + max_large_obj = abs_cost; + } + } + if (num_large_obj > 0) + { + fprintf(stderr, + "WARNING: %d objective %s |value| >= %.1e (largest %.3e); the problem is badly scaled.\n", + num_large_obj, + (num_large_obj == 1 ? "coefficient has" : "coefficients have"), + OBJECTIVE_LARGE_VALUE, + max_large_obj); + } + filter_constraint_matrix_entries(&working, original, params); return working; } @@ -434,6 +513,16 @@ void free_preprocessed_problem(const lp_problem_t *preprocessed, const lp_proble { if (preprocessed->objective_vector != original->objective_vector) free(preprocessed->objective_vector); + if (preprocessed->variable_lower_bound != original->variable_lower_bound) + { + free(preprocessed->variable_lower_bound); + free(preprocessed->variable_upper_bound); + } + if (preprocessed->constraint_lower_bound != original->constraint_lower_bound) + { + free(preprocessed->constraint_lower_bound); + free(preprocessed->constraint_upper_bound); + } if (preprocessed->constraint_matrix_values != original->constraint_matrix_values) { free(preprocessed->constraint_matrix_row_pointers); @@ -536,6 +625,7 @@ void print_initial_info(const pdhg_parameters_t *params, const lp_problem_t *pro default_params.termination_criteria.eps_feas_polish_relative); PRINT_DIFF_BOOL("presolve", params->presolve, default_params.presolve); PRINT_DIFF_DBL("matrix_zero_tol", params->matrix_zero_tol, default_params.matrix_zero_tol); + PRINT_DIFF_DBL("infinite_bound", params->infinite_bound, default_params.infinite_bound); } #undef PRINT_DIFF_INT