Skip to content

Commit 3e52ffe

Browse files
nguidottinguidotti
committed
fix infeasible list (NVIDIA#694)
1 parent 426b445 commit 3e52ffe

File tree

7 files changed

+118
-37
lines changed

7 files changed

+118
-37
lines changed

‎cpp/src/dual_simplex/basis_solves.cpp‎

Lines changed: 14 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -613,6 +613,8 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
613613
template <typename i_t, typename f_t>
614614
i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
615615
const simplex_solver_settings_t<i_t, f_t>& settings,
616+
const std::vector<f_t>& lower,
617+
const std::vector<f_t>& upper,
616618
const std::vector<i_t>& deficient,
617619
const std::vector<i_t>& slacks_needed,
618620
std::vector<i_t>& basis_list,
@@ -658,7 +660,15 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
658660
nonbasic_list[nonbasic_map[replace_j]] = bad_j;
659661
vstatus[replace_j] = variable_status_t::BASIC;
660662
// This is the main issue. What value should bad_j take on.
661-
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
663+
if (lower[bad_j] == -inf && upper[bad_j] == inf) {
664+
vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
665+
} else if (lower[bad_j] > -inf) {
666+
vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
667+
} else if (upper[bad_j] < inf) {
668+
vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
669+
} else {
670+
assert(1 == 0);
671+
}
662672
}
663673

664674
return 0;
@@ -849,6 +859,8 @@ template int factorize_basis<int>(const csc_matrix_t<int, double>& A,
849859

850860
template int basis_repair<int, double>(const csc_matrix_t<int, double>& A,
851861
const simplex_solver_settings_t<int, double>& settings,
862+
const std::vector<double>& lower,
863+
const std::vector<double>& upper,
852864
const std::vector<int>& deficient,
853865
const std::vector<int>& slacks_needed,
854866
std::vector<int>& basis_list,

‎cpp/src/dual_simplex/basis_solves.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -42,6 +42,8 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
4242
template <typename i_t, typename f_t>
4343
i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
4444
const simplex_solver_settings_t<i_t, f_t>& settings,
45+
const std::vector<f_t>& lower,
46+
const std::vector<f_t>& upper,
4547
const std::vector<i_t>& deficient,
4648
const std::vector<i_t>& slacks_needed,
4749
std::vector<i_t>& basis_list,

‎cpp/src/dual_simplex/basis_updates.cpp‎

Lines changed: 5 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -2046,6 +2046,8 @@ template <typename i_t, typename f_t>
20462046
int basis_update_mpf_t<i_t, f_t>::refactor_basis(
20472047
const csc_matrix_t<i_t, f_t>& A,
20482048
const simplex_solver_settings_t<i_t, f_t>& settings,
2049+
const std::vector<f_t>& lower,
2050+
const std::vector<f_t>& upper,
20492051
std::vector<i_t>& basic_list,
20502052
std::vector<i_t>& nonbasic_list,
20512053
std::vector<variable_status_t>& vstatus)
@@ -2066,7 +2068,8 @@ int basis_update_mpf_t<i_t, f_t>::refactor_basis(
20662068
deficient,
20672069
slacks_needed) == -1) {
20682070
settings.log.debug("Initial factorization failed\n");
2069-
basis_repair(A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
2071+
basis_repair(
2072+
A, settings, lower, upper, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
20702073

20712074
#ifdef CHECK_BASIS_REPAIR
20722075
const i_t m = A.m;

‎cpp/src/dual_simplex/basis_updates.hpp‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -373,6 +373,8 @@ class basis_update_mpf_t {
373373
// Compute L*U = A(p, basic_list)
374374
int refactor_basis(const csc_matrix_t<i_t, f_t>& A,
375375
const simplex_solver_settings_t<i_t, f_t>& settings,
376+
const std::vector<f_t>& lower,
377+
const std::vector<f_t>& upper,
376378
std::vector<i_t>& basic_list,
377379
std::vector<i_t>& nonbasic_list,
378380
std::vector<variable_status_t>& vstatus);

‎cpp/src/dual_simplex/crossover.cpp‎

Lines changed: 28 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -785,8 +785,15 @@ i_t primal_push(const lp_problem_t<i_t, f_t>& lp,
785785
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
786786
if (rank != m) {
787787
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
788-
basis_repair(
789-
lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
788+
basis_repair(lp.A,
789+
settings,
790+
lp.lower,
791+
lp.upper,
792+
deficient,
793+
slacks_needed,
794+
basic_list,
795+
nonbasic_list,
796+
vstatus);
790797
if (factorize_basis(
791798
lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
792799
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
@@ -1132,7 +1139,15 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
11321139
rank = factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
11331140
if (rank != m) {
11341141
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
1135-
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
1142+
basis_repair(lp.A,
1143+
settings,
1144+
lp.lower,
1145+
lp.upper,
1146+
deficient,
1147+
slacks_needed,
1148+
basic_list,
1149+
nonbasic_list,
1150+
vstatus);
11361151
if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
11371152
-1) {
11381153
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);
@@ -1323,7 +1338,15 @@ crossover_status_t crossover(const lp_problem_t<i_t, f_t>& lp,
13231338
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
13241339
if (rank != m) {
13251340
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
1326-
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
1341+
basis_repair(lp.A,
1342+
settings,
1343+
lp.lower,
1344+
lp.upper,
1345+
deficient,
1346+
slacks_needed,
1347+
basic_list,
1348+
nonbasic_list,
1349+
vstatus);
13271350
if (factorize_basis(
13281351
lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
13291352
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);

‎cpp/src/dual_simplex/phase2.cpp‎

Lines changed: 55 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -623,14 +623,17 @@ f_t compute_initial_primal_infeasibilities(const lp_problem_t<i_t, f_t>& lp,
623623
const std::vector<i_t>& basic_list,
624624
const std::vector<f_t>& x,
625625
std::vector<f_t>& squared_infeasibilities,
626-
std::vector<i_t>& infeasibility_indices)
626+
std::vector<i_t>& infeasibility_indices,
627+
f_t& primal_inf)
627628
{
628629
const i_t m = lp.num_rows;
629630
const i_t n = lp.num_cols;
630-
squared_infeasibilities.resize(n, 0.0);
631+
squared_infeasibilities.resize(n);
632+
std::fill(squared_infeasibilities.begin(), squared_infeasibilities.end(), 0.0);
631633
infeasibility_indices.reserve(n);
632634
infeasibility_indices.clear();
633-
f_t primal_inf = 0.0;
635+
f_t primal_inf_squared = 0.0;
636+
primal_inf = 0.0;
634637
for (i_t k = 0; k < m; ++k) {
635638
const i_t j = basic_list[k];
636639
const f_t lower_infeas = lp.lower[j] - x[j];
@@ -640,10 +643,11 @@ f_t compute_initial_primal_infeasibilities(const lp_problem_t<i_t, f_t>& lp,
640643
const f_t square_infeas = infeas * infeas;
641644
squared_infeasibilities[j] = square_infeas;
642645
infeasibility_indices.push_back(j);
643-
primal_inf += square_infeas;
646+
primal_inf_squared += square_infeas;
647+
primal_inf += infeas;
644648
}
645649
}
646-
return primal_inf;
650+
return primal_inf_squared;
647651
}
648652

649653
template <typename i_t, typename f_t>
@@ -2241,7 +2245,8 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
22412245
assert(superbasic_list.size() == 0);
22422246
assert(nonbasic_list.size() == n - m);
22432247

2244-
if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
2248+
if (ft.refactor_basis(lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) >
2249+
0) {
22452250
return dual::status_t::NUMERICAL;
22462251
}
22472252

@@ -2268,7 +2273,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
22682273

22692274
#ifdef COMPUTE_DUAL_RESIDUAL
22702275
std::vector<f_t> dual_res1;
2271-
compute_dual_residual(lp.A, objective, y, z, dual_res1);
2276+
phase2::compute_dual_residual(lp.A, objective, y, z, dual_res1);
22722277
f_t dual_res_norm = vector_norm_inf<i_t, f_t>(dual_res1);
22732278
if (dual_res_norm > settings.tight_tol) {
22742279
settings.log.printf("|| A'*y + z - c || %e\n", dual_res_norm);
@@ -2357,8 +2362,15 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
23572362
std::vector<uint8_t> bounded_variables(n, 0);
23582363
phase2::compute_bounded_info(lp.lower, lp.upper, bounded_variables);
23592364

2360-
f_t primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
2361-
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
2365+
f_t primal_infeasibility;
2366+
f_t primal_infeasibility_squared =
2367+
phase2::compute_initial_primal_infeasibilities(lp,
2368+
settings,
2369+
basic_list,
2370+
x,
2371+
squared_infeasibilities,
2372+
infeasibility_indices,
2373+
primal_infeasibility);
23622374

23632375
#ifdef CHECK_BASIC_INFEASIBILITIES
23642376
phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 0);
@@ -2556,9 +2568,15 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
25562568
std::vector<f_t> unperturbed_x(n);
25572569
phase2::compute_primal_solution_from_basis(
25582570
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
2559-
x = unperturbed_x;
2560-
primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
2561-
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
2571+
x = unperturbed_x;
2572+
primal_infeasibility_squared =
2573+
phase2::compute_initial_primal_infeasibilities(lp,
2574+
settings,
2575+
basic_list,
2576+
x,
2577+
squared_infeasibilities,
2578+
infeasibility_indices,
2579+
primal_infeasibility);
25622580
settings.log.printf("Updated primal infeasibility: %e\n", primal_infeasibility);
25632581

25642582
objective = lp.objective;
@@ -2593,9 +2611,15 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
25932611
std::vector<f_t> unperturbed_x(n);
25942612
phase2::compute_primal_solution_from_basis(
25952613
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
2596-
x = unperturbed_x;
2597-
primal_infeasibility = phase2::compute_initial_primal_infeasibilities(
2598-
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
2614+
x = unperturbed_x;
2615+
primal_infeasibility_squared =
2616+
phase2::compute_initial_primal_infeasibilities(lp,
2617+
settings,
2618+
basic_list,
2619+
x,
2620+
squared_infeasibilities,
2621+
infeasibility_indices,
2622+
primal_infeasibility);
25992623

26002624
const f_t orig_dual_infeas = phase2::dual_infeasibility(
26012625
lp, settings, vstatus, z, settings.tight_tol, settings.dual_tol);
@@ -2810,7 +2834,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
28102834
delta_xB_0_sparse.i,
28112835
squared_infeasibilities,
28122836
infeasibility_indices,
2813-
primal_infeasibility);
2837+
primal_infeasibility_squared);
28142838
// Update primal infeasibilities due to changes in basic variables
28152839
// from the leaving and entering variables
28162840
phase2::update_primal_infeasibilities(lp,
@@ -2822,7 +2846,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
28222846
scaled_delta_xB_sparse.i,
28232847
squared_infeasibilities,
28242848
infeasibility_indices,
2825-
primal_infeasibility);
2849+
primal_infeasibility_squared);
28262850
// Update the entering variable
28272851
phase2::update_single_primal_infeasibility(lp.lower,
28282852
lp.upper,
@@ -2883,14 +2907,15 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
28832907
#endif
28842908
if (should_refactor) {
28852909
bool should_recompute_x = false;
2886-
if (ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus) > 0) {
2910+
if (ft.refactor_basis(
2911+
lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus) > 0) {
28872912
should_recompute_x = true;
28882913
settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
28892914
if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
28902915
i_t count = 0;
28912916
i_t deficient_size;
2892-
while ((deficient_size =
2893-
ft.refactor_basis(lp.A, settings, basic_list, nonbasic_list, vstatus)) > 0) {
2917+
while ((deficient_size = ft.refactor_basis(
2918+
lp.A, settings, lp.lower, lp.upper, basic_list, nonbasic_list, vstatus)) > 0) {
28942919
settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
28952920
iter,
28962921
static_cast<int>(deficient_size));
@@ -2912,8 +2937,14 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
29122937
lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
29132938
x = unperturbed_x;
29142939
}
2915-
phase2::compute_initial_primal_infeasibilities(
2916-
lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
2940+
primal_infeasibility_squared =
2941+
phase2::compute_initial_primal_infeasibilities(lp,
2942+
settings,
2943+
basic_list,
2944+
x,
2945+
squared_infeasibilities,
2946+
infeasibility_indices,
2947+
primal_infeasibility);
29172948
}
29182949
#ifdef CHECK_BASIC_INFEASIBILITIES
29192950
phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 7);
@@ -2951,7 +2982,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase,
29512982
iter,
29522983
compute_user_objective(lp, obj),
29532984
infeasibility_indices.size(),
2954-
primal_infeasibility,
2985+
primal_infeasibility_squared,
29552986
sum_perturb,
29562987
now);
29572988
}

‎cpp/src/dual_simplex/primal.cpp‎

Lines changed: 10 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
/* clang-format off */
22
/*
3-
* SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
3+
* SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
44
* SPDX-License-Identifier: Apache-2.0
55
*/
66
/* clang-format on */
@@ -298,7 +298,15 @@ primal::status_t primal_phase2(i_t phase,
298298
factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed);
299299
if (rank != m) {
300300
settings.log.debug("Failed to factorize basis. rank %d m %d\n", rank, m);
301-
basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
301+
basis_repair(lp.A,
302+
settings,
303+
lp.lower,
304+
lp.upper,
305+
deficient,
306+
slacks_needed,
307+
basic_list,
308+
nonbasic_list,
309+
vstatus);
302310
if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
303311
-1) {
304312
settings.log.printf("Failed to factorize basis after repair. rank %d m %d\n", rank, m);

0 commit comments

Comments
 (0)