diff --git a/cpp/src/barrier/barrier.cu b/cpp/src/barrier/barrier.cu index bf7ffbdf86..a9b3aa891a 100644 --- a/cpp/src/barrier/barrier.cu +++ b/cpp/src/barrier/barrier.cu @@ -3000,14 +3000,14 @@ i_t barrier_solver_t::gpu_compute_search_direction(iteration_data_t %e (solve_err=%e)\n", old_dp, dual_perturb, solve_err); } else if (solve_err < 1e-4) { f_t old_dp = dual_perturb; dual_perturb = std::max(min_perturb, dual_perturb / 10.0); primal_perturb = std::max(min_perturb, primal_perturb / 10.0); if (old_dp != dual_perturb) { - settings.log.printf( + settings.log.debug( " reg DOWN: %e -> %e (solve_err=%e)\n", old_dp, dual_perturb, solve_err); } } diff --git a/cpp/src/barrier/translate_soc.hpp b/cpp/src/barrier/translate_soc.hpp index 451209bed8..3209dee3c6 100644 --- a/cpp/src/barrier/translate_soc.hpp +++ b/cpp/src/barrier/translate_soc.hpp @@ -10,8 +10,10 @@ #include #include +#include #include #include +#include #include #include @@ -64,6 +66,27 @@ void convert_quadratic_constraints_to_second_order_cones( // Use a practical tolerance for text-parsed MPS numeric values. const f_t tol = std::numeric_limits::epsilon() * 2; + // Derive implied lower bounds from singleton inequality rows. + // Used to check if SOC head variables have implied non-negativity from the constraint system + // without actually modifying the variable bounds (which would add barrier terms). + std::vector implied_lower(n, -std::numeric_limits::infinity()); + for (i_t i = 0; i < csr_A.m; i++) { + const i_t row_start = csr_A.row_start[i]; + const i_t row_end = csr_A.row_start[i + 1]; + if (row_end - row_start != 1) { continue; } + const i_t j = csr_A.j[row_start]; + const f_t a = csr_A.x[row_start]; + const f_t b = user_problem.rhs[i]; + const char sense = user_problem.row_sense[i]; + if (std::abs(a) < tol) { continue; } + const f_t bound = b / a; + if (sense == 'G' && a > 0) { + implied_lower[j] = std::max(implied_lower[j], bound); + } else if (sense == 'L' && a < 0) { + implied_lower[j] = std::max(implied_lower[j], bound); + } + } + // SOC conversion accepts: // 1) diagonal Lorentz-form QCMATRIX rows: // -s*x_head^2 + sum_i s*x_tail_i^2 <= 0 (any common s > 0; divide by s to normalize) @@ -121,9 +144,6 @@ void convert_quadratic_constraints_to_second_order_cones( "Quadratic constraint '%s' ROWS type must be 'L' (<=) or 'G' (>=)", qc.constraint_row_name.c_str()); normalize_quadratic_constraint_greater_to_less(qc); - cuopt_expects((qc.rhs_value < tol) && (qc.rhs_value > -tol), - error_type_t::ValidationError, - "SOC conversion currently requires rhs = 0 for quadratic constraints"); cuopt_expects(qc.linear_values.size() == qc.linear_indices.size(), error_type_t::ValidationError, "Quadratic constraint '%s' linear_values and linear_indices length mismatch", @@ -178,7 +198,7 @@ void convert_quadratic_constraints_to_second_order_cones( return std::abs(a - b) <= tol * scale; }; - // Sort COO by (row, col); O(nnz log nnz). Enforce at most one stored entry per row (SOC CSR). + // Sort COO by (row, col); O(nnz log nnz). std::vector perm(q_nnz); std::iota(perm.begin(), perm.end(), size_t{0}); std::sort(perm.begin(), perm.end(), [&](size_t a, size_t b) { @@ -190,6 +210,7 @@ void convert_quadratic_constraints_to_second_order_cones( std::vector> q_entries; q_entries.reserve(q_nnz); + bool has_duplicate_rows = false; for (size_t t = 0; t < static_cast(q_nnz); ++t) { const size_t ix = perm[t]; const i_t r = qc.rows[ix]; @@ -204,12 +225,7 @@ void convert_quadratic_constraints_to_second_order_cones( static_cast(n)); if (!q_entries.empty()) { const i_t prev_r = std::get<0>(q_entries.back()); - cuopt_expects(r != prev_r, - error_type_t::ValidationError, - "Quadratic constraint '%s' Q row %d: expected at most one stored entry per " - "row (CSR layout); duplicate or unsorted row in COO", - qc.constraint_row_name.c_str(), - static_cast(r)); + if (r == prev_r) { has_duplicate_rows = true; } } q_entries.emplace_back(r, c, v); } @@ -221,6 +237,7 @@ void convert_quadratic_constraints_to_second_order_cones( neg_diag_rows.reserve(1); offdiag_entries.reserve(4); + bool has_near_zero_diag = false; for (const auto& [r, c, v] : q_entries) { if (r == c) { if (v > tol) { @@ -228,13 +245,7 @@ void convert_quadratic_constraints_to_second_order_cones( } else if (v < -tol) { neg_diag_rows.emplace_back(r, v); } else { - cuopt_expects(false, - error_type_t::ValidationError, - "Quadratic constraint '%s' Q row %d: diagonal SOC entry is near zero " - "(%.17g)", - qc.constraint_row_name.c_str(), - static_cast(r), - static_cast(v)); + has_near_zero_diag = true; } } else { offdiag_entries.emplace_back(r, c, v); @@ -247,6 +258,28 @@ void convert_quadratic_constraints_to_second_order_cones( tail_vars.push_back(pr.first); } + // Determine whether to use the general convex quadratic path. + // The general path is needed when Q does not fit any special SOC pattern, + // or when the RHS is nonzero (special cases require rhs = 0). + const bool has_nonzero_rhs = !(qc.rhs_value < tol && qc.rhs_value > -tol); + bool has_nonuniform_diag = false; + if (pos_diag_rows.size() > 1) { + const f_t first_val = pos_diag_rows[0].second; + for (size_t k = 1; k < pos_diag_rows.size(); k++) { + const f_t scale = + std::max({f_t(1), std::abs(first_val), std::abs(pos_diag_rows[k].second)}); + if (std::abs(pos_diag_rows[k].second - first_val) > tol * scale) { + has_nonuniform_diag = true; + break; + } + } + } + const bool use_general_path = has_duplicate_rows || has_near_zero_diag || has_nonzero_rhs || + has_nonuniform_diag || offdiag_entries.size() > 2 || + (offdiag_entries.size() == 1) || (neg_diag_rows.size() > 1) || + (!neg_diag_rows.empty() && has_linear_part) || + (!neg_diag_rows.empty() && !offdiag_entries.empty()); + f_t uniform_s = 0; bool have_uniform_s = false; auto note_positive_s = [&](f_t v) { @@ -276,228 +309,504 @@ void convert_quadratic_constraints_to_second_order_cones( char is_rotated = 0; i_t head = -1; - if (offdiag_entries.empty()) { - if (!has_linear_part) { - if (pos_diag_rows.empty()) { - cuopt_expects(neg_diag_rows.size() == 1 && q_nnz == 1, - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: expected tail diagonals +s with head -s, " - "or a single head row with q_nnz=1", - qc.constraint_row_name.c_str()); - const f_t neg_v = neg_diag_rows[0].second; - cuopt_expects(neg_v < -tol, - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: cone head diagonal must be negative " - "(%.17g)", - qc.constraint_row_name.c_str(), - static_cast(neg_v)); - uniform_s = -neg_v; - have_uniform_s = true; - head = neg_diag_rows[0].first; + if (!use_general_path) { + // Special-case rhs == 0 requirement for SOC patterns + cuopt_expects((qc.rhs_value < tol) && (qc.rhs_value > -tol), + error_type_t::ValidationError, + "SOC conversion currently requires rhs = 0 for quadratic constraints " + "(constraint '%s' has rhs %.17g)", + qc.constraint_row_name.c_str(), + static_cast(qc.rhs_value)); + + if (offdiag_entries.empty()) { + if (!has_linear_part) { + if (pos_diag_rows.empty()) { + cuopt_expects( + neg_diag_rows.size() == 1 && q_nnz == 1, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: expected tail diagonals +s with head -s, " + "or a single head row with q_nnz=1", + qc.constraint_row_name.c_str()); + const f_t neg_v = neg_diag_rows[0].second; + cuopt_expects(neg_v < -tol, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: cone head diagonal must be negative " + "(%.17g)", + qc.constraint_row_name.c_str(), + static_cast(neg_v)); + uniform_s = -neg_v; + have_uniform_s = true; + head = neg_diag_rows[0].first; + cuopt_expects( + static_cast(tail_vars.size()) == q_nnz - 1, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: expected %d diagonal +s entries (tails), found %zu", + qc.constraint_row_name.c_str(), + static_cast(q_nnz - 1), + tail_vars.size()); + cone.reserve(1); + cone.push_back(head); + cone_dim = static_cast(cone.size()); + is_rotated = 0; + } else { + for (const std::pair& pr : pos_diag_rows) { + note_positive_s(pr.second); + } + cuopt_expects( + have_uniform_s, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: could not infer uniform positive scale s", + qc.constraint_row_name.c_str()); + cuopt_expects( + neg_diag_rows.size() == 1, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: expected exactly one diagonal -s (cone head) for " + "%zu tail entries, found %zu negative diagonals", + qc.constraint_row_name.c_str(), + tail_vars.size(), + neg_diag_rows.size()); + cuopt_expects( + static_cast(tail_vars.size()) == q_nnz - 1, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: expected %d diagonal +s entries (tails), found %zu", + qc.constraint_row_name.c_str(), + static_cast(q_nnz - 1), + tail_vars.size()); + const f_t neg_v = neg_diag_rows[0].second; + cuopt_expects( + approx_eq_scaled(neg_v, -uniform_s), + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC Q: cone head diagonal must be -s with the same s as " + "positive tail diagonals; head %.17g vs -s = %.17g", + qc.constraint_row_name.c_str(), + static_cast(neg_v), + static_cast(-uniform_s)); + head = neg_diag_rows[0].first; + // The SOC ||tail|| <= head requires head >= 0. Check explicit bound + // or implied bound from singleton inequality constraints. + cuopt_expects(std::max(user_problem.lower[head], implied_lower[head]) >= 0, + error_type_t::ValidationError, + "Quadratic constraint '%s' SOC head variable (index %d) must have a " + "non-negative lower bound for the constraint to be convex", + qc.constraint_row_name.c_str(), + static_cast(head)); + cone.reserve(q_nnz); + cone.push_back(head); + cone.insert(cone.end(), tail_vars.begin(), tail_vars.end()); + cone_dim = static_cast(cone.size()); + is_rotated = 0; + } + } else { cuopt_expects( - static_cast(tail_vars.size()) == q_nnz - 1, + neg_diag_rows.empty(), error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: expected %d diagonal +s entries (tails), found %zu", - qc.constraint_row_name.c_str(), - static_cast(q_nnz - 1), - tail_vars.size()); - cone.reserve(1); - cone.push_back(head); - cone_dim = static_cast(cone.size()); - is_rotated = 0; - } else { + "Quadratic constraint '%s' with linear terms cannot contain negative diagonal " + "Q entries", + qc.constraint_row_name.c_str()); + cuopt_expects(affine_head >= 0, + error_type_t::ValidationError, + "Quadratic constraint '%s' internal error: affine SOC head index invalid", + qc.constraint_row_name.c_str()); for (const std::pair& pr : pos_diag_rows) { note_positive_s(pr.second); } cuopt_expects(have_uniform_s, error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: could not infer uniform positive scale s", + "Quadratic constraint '%s' with linear terms must have at least one " + "diagonal +s term in Q", qc.constraint_row_name.c_str()); - cuopt_expects( - neg_diag_rows.size() == 1, - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: expected exactly one diagonal -s (cone head) for " - "%zu tail entries, found %zu negative diagonals", - qc.constraint_row_name.c_str(), - tail_vars.size(), - neg_diag_rows.size()); - cuopt_expects( - static_cast(tail_vars.size()) == q_nnz - 1, - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: expected %d diagonal +s entries (tails), found %zu", - qc.constraint_row_name.c_str(), - static_cast(q_nnz - 1), - tail_vars.size()); - const f_t neg_v = neg_diag_rows[0].second; - cuopt_expects( - approx_eq_scaled(neg_v, -uniform_s), - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: cone head diagonal must be -s with the same s as " - "positive tail diagonals; head %.17g vs -s = %.17g", - qc.constraint_row_name.c_str(), - static_cast(neg_v), - static_cast(-uniform_s)); - head = neg_diag_rows[0].first; - cone.reserve(q_nnz); - cone.push_back(head); + cuopt_expects(!tail_vars.empty(), + error_type_t::ValidationError, + "Quadratic constraint '%s' with linear terms must have at least one " + "diagonal +s term in Q", + qc.constraint_row_name.c_str()); + for (const i_t tail : tail_vars) { + cuopt_expects( + tail != affine_head, + error_type_t::ValidationError, + "Quadratic constraint '%s' with linear terms requires the linear head variable to be " + "distinct from quadratic diagonal variables", + qc.constraint_row_name.c_str()); + } + + cone.reserve(tail_vars.size() + 1); + cone.push_back(affine_head); cone.insert(cone.end(), tail_vars.begin(), tail_vars.end()); - cone_dim = static_cast(cone.size()); - is_rotated = 0; + cone_dim = static_cast(tail_vars.size() + 2); + is_rotated = 1; + rotated_cones.push_back(rotated_soc_t{affine_head, -1, tail_vars, true, 1}); } } else { - cuopt_expects( - neg_diag_rows.empty(), - error_type_t::ValidationError, - "Quadratic constraint '%s' with linear terms cannot contain negative diagonal " - "Q entries", - qc.constraint_row_name.c_str()); - cuopt_expects(affine_head >= 0, + cuopt_expects(!has_linear_part, error_type_t::ValidationError, - "Quadratic constraint '%s' internal error: affine SOC head index invalid", + "Quadratic constraint '%s' with linear terms cannot include rotated-SOC " + "off-diagonal entries", qc.constraint_row_name.c_str()); + cuopt_expects(neg_diag_rows.empty(), + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q cannot contain diagonal head " + "entries; found %zu negative diagonals", + qc.constraint_row_name.c_str(), + neg_diag_rows.size()); for (const std::pair& pr : pos_diag_rows) { note_positive_s(pr.second); } cuopt_expects(have_uniform_s, error_type_t::ValidationError, - "Quadratic constraint '%s' with linear terms must have at least one " - "diagonal +s term in Q", + "Quadratic constraint '%s' rotated SOC Q: could not infer uniform scale s", + qc.constraint_row_name.c_str()); + cuopt_expects( + offdiag_entries.size() == 2, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q must contain exactly one symmetric off-diagonal " + "pair (-d,-d); found %zu off-diagonal entries", + qc.constraint_row_name.c_str(), + offdiag_entries.size()); + + const i_t a = std::get<0>(offdiag_entries[0]); + const i_t b = std::get<1>(offdiag_entries[0]); + const f_t v0 = std::get<2>(offdiag_entries[0]); + cuopt_expects( + v0 < -tol, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g", + qc.constraint_row_name.c_str(), + static_cast(v0)); + cuopt_expects(a != b, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q off-diagonal pair must use distinct " + "variables", qc.constraint_row_name.c_str()); - cuopt_expects(!tail_vars.empty(), + cuopt_expects(std::get<0>(offdiag_entries[1]) == b && std::get<1>(offdiag_entries[1]) == a, error_type_t::ValidationError, - "Quadratic constraint '%s' with linear terms must have at least one " - "diagonal +s term in Q", + "Quadratic constraint '%s' rotated SOC Q must have symmetric entries (a,b) " + "and (b,a) with the same value", + qc.constraint_row_name.c_str()); + const f_t v1 = std::get<2>(offdiag_entries[1]); + cuopt_expects( + v1 < -tol, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g", + qc.constraint_row_name.c_str(), + static_cast(v1)); + cuopt_expects( + approx_eq_scaled(v0, v1), + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q symmetric off-diagonals must match; got %.17g " + "and %.17g", + qc.constraint_row_name.c_str(), + static_cast(v0), + static_cast(v1)); + const f_t cross_d = -v0; + cuopt_expects( + cross_d > tol, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q cross coefficient d = -Q_off must be positive", + qc.constraint_row_name.c_str()); + const f_t head_lift_sqrt_ratio = std::sqrt(cross_d / uniform_s); + cuopt_expects(std::isfinite(static_cast(head_lift_sqrt_ratio)), + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q head lift ratio sqrt(d/s) is not " + "finite (d=%.17g, s=%.17g)", + qc.constraint_row_name.c_str(), + static_cast(cross_d), + static_cast(uniform_s)); + cuopt_expects(static_cast(tail_vars.size()) == q_nnz - 2, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q: expected %d diagonal +s entries " + "(tails), found %zu", + qc.constraint_row_name.c_str(), + static_cast(q_nnz - 2), + tail_vars.size()); + cuopt_expects(q_nnz >= 3, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC Q must have at least 1 tail entry", qc.constraint_row_name.c_str()); - for (const i_t tail : tail_vars) { - cuopt_expects( - tail != affine_head, - error_type_t::ValidationError, - "Quadratic constraint '%s' with linear terms requires the linear head variable to be " - "distinct from quadratic diagonal variables", - qc.constraint_row_name.c_str()); - } - cone.reserve(tail_vars.size() + 1); - cone.push_back(affine_head); + cone.reserve(q_nnz); + cone.push_back(a); + cone.push_back(b); cone.insert(cone.end(), tail_vars.begin(), tail_vars.end()); - cone_dim = static_cast(tail_vars.size() + 2); + cone_dim = static_cast(cone.size()); is_rotated = 1; - rotated_cones.push_back(rotated_soc_t{affine_head, -1, tail_vars, true, 1}); - } - } else { - cuopt_expects(!has_linear_part, - error_type_t::ValidationError, - "Quadratic constraint '%s' with linear terms cannot include rotated-SOC " - "off-diagonal entries", - qc.constraint_row_name.c_str()); - cuopt_expects(neg_diag_rows.empty(), - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q cannot contain diagonal head " - "entries; found %zu negative diagonals", - qc.constraint_row_name.c_str(), - neg_diag_rows.size()); - for (const std::pair& pr : pos_diag_rows) { - note_positive_s(pr.second); + // Rotated SOC ||tail||^2 <= 2*a*b requires a >= 0 and b >= 0. + // Check explicit bound or implied bound from singleton inequality constraints. + cuopt_expects(std::max(user_problem.lower[a], implied_lower[a]) >= 0, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC head variable (index %d) must have a " + "non-negative lower bound for the constraint to be convex", + qc.constraint_row_name.c_str(), + static_cast(a)); + cuopt_expects(std::max(user_problem.lower[b], implied_lower[b]) >= 0, + error_type_t::ValidationError, + "Quadratic constraint '%s' rotated SOC head variable (index %d) must have a " + "non-negative lower bound for the constraint to be convex", + qc.constraint_row_name.c_str(), + static_cast(b)); + rotated_cones.push_back(rotated_soc_t{a, b, tail_vars, false, head_lift_sqrt_ratio}); } - cuopt_expects(have_uniform_s, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q: could not infer uniform scale s", - qc.constraint_row_name.c_str()); - cuopt_expects( - offdiag_entries.size() == 2, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q must contain exactly one symmetric off-diagonal " - "pair (-d,-d); found %zu off-diagonal entries", - qc.constraint_row_name.c_str(), - offdiag_entries.size()); - - const i_t a = std::get<0>(offdiag_entries[0]); - const i_t b = std::get<1>(offdiag_entries[0]); - const f_t v0 = std::get<2>(offdiag_entries[0]); - cuopt_expects( - v0 < -tol, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g", - qc.constraint_row_name.c_str(), - static_cast(v0)); - cuopt_expects(a != b, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q off-diagonal pair must use distinct " - "variables", - qc.constraint_row_name.c_str()); - cuopt_expects(std::get<0>(offdiag_entries[1]) == b && std::get<1>(offdiag_entries[1]) == a, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q must have symmetric entries (a,b) " - "and (b,a) with the same value", - qc.constraint_row_name.c_str()); - const f_t v1 = std::get<2>(offdiag_entries[1]); - cuopt_expects( - v1 < -tol, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q off-diagonal must be negative; got %.17g", - qc.constraint_row_name.c_str(), - static_cast(v1)); - cuopt_expects( - approx_eq_scaled(v0, v1), - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q symmetric off-diagonals must match; got %.17g " - "and %.17g", - qc.constraint_row_name.c_str(), - static_cast(v0), - static_cast(v1)); - const f_t cross_d = -v0; - cuopt_expects( - cross_d > tol, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q cross coefficient d = -Q_off must be positive", - qc.constraint_row_name.c_str()); - const f_t head_lift_sqrt_ratio = std::sqrt(cross_d / uniform_s); - cuopt_expects(std::isfinite(static_cast(head_lift_sqrt_ratio)), + + cuopt_expects(have_uniform_s && uniform_s > tol, error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q head lift ratio sqrt(d/s) is not " - "finite (d=%.17g, s=%.17g)", + "Quadratic constraint '%s' SOC Q: uniform scale s must be positive (got %.17g)", qc.constraint_row_name.c_str(), - static_cast(cross_d), static_cast(uniform_s)); - cuopt_expects(static_cast(tail_vars.size()) == q_nnz - 2, - error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q: expected %d diagonal +s entries " - "(tails), found %zu", - qc.constraint_row_name.c_str(), - static_cast(q_nnz - 2), - tail_vars.size()); - cuopt_expects(q_nnz >= 3, + qc_soc_uniform_scale[qc_i] = uniform_s; + + for (const i_t var : cone) { + cuopt_expects(var >= 0 && var < static_cast(is_cone_var.size()), + error_type_t::ValidationError, + "SOC variable index %d is outside [0, %zu)", + static_cast(var), + is_cone_var.size()); + } + cone_dims.push_back(cone_dim); + cone_vars.push_back(std::move(cone)); + cone_is_rotated.push_back(is_rotated); + + } else { + // ========================================================================= + // General convex quadratic constraint path: + // x^T Q x + c^T x <= alpha + + // Invalidate affine head for this QC — the general path handles the linear part directly + qc_affine_heads[qc_i] = -1; + // where Q is (possibly unsymmetric) and H = Q + Q^T must be PSD. + // ========================================================================= + const f_t alpha = qc.rhs_value; + + // Step 1: Build H such that (1/2) x^T H x equals the quadratic form sum_k + // v_k*x_{r_k}*x_{c_k}. For diagonal entry (r,r,v): H(r,r) += 2*v (since (1/2)*H(r,r)*x_r^2 = + // v*x_r^2) For off-diagonal entry (r,c,v): H(max,min) += v (since + // (1/2)*(H(r,c)+H(c,r))*x_r*x_c = v*x_r*x_c) Store lower triangle only in CSC. + // + // Use a dense accumulator indexed by the variables appearing in Q. + + // Collect distinct variable indices and build local-to-global mapping + std::vector var_set; + var_set.reserve(2 * q_nnz); + std::vector global_to_local(n, -1); + for (size_t t = 0; t < static_cast(q_nnz); ++t) { + const i_t r = qc.rows[t]; + const i_t c = qc.cols[t]; + if (global_to_local[r] == -1) { + global_to_local[r] = static_cast(var_set.size()); + var_set.push_back(r); + } + if (global_to_local[c] == -1) { + global_to_local[c] = static_cast(var_set.size()); + var_set.push_back(c); + } + } + const i_t n_local = static_cast(var_set.size()); + + // Dense lower-triangle accumulator (column-major: H_dense[col * n_local + row] for row >= + // col) + std::vector H_dense(n_local * n_local, f_t(0)); + for (size_t t = 0; t < static_cast(q_nnz); ++t) { + const i_t r = global_to_local[qc.rows[t]]; + const i_t c = global_to_local[qc.cols[t]]; + const f_t v = qc.vals[t]; + if (r == c) { + H_dense[c * n_local + r] += f_t(2) * v; + } else { + const i_t hi = std::max(r, c); + const i_t hj = std::min(r, c); + H_dense[hj * n_local + hi] += v; + } + } + + // Gather nonzeros from dense accumulator into CSC (lower triangle, local indices) + i_t h_nnz = 0; + for (i_t j = 0; j < n_local; j++) { + for (i_t i = j; i < n_local; i++) { + if (H_dense[j * n_local + i] != f_t(0)) { h_nnz++; } + } + } + + dual_simplex::csc_matrix_t H_csc(n_local, n_local, h_nnz); + { + i_t p = 0; + for (i_t j = 0; j < n_local; j++) { + H_csc.col_start[j] = p; + for (i_t i = j; i < n_local; i++) { + const f_t val = H_dense[j * n_local + i]; + if (val != f_t(0)) { + H_csc.i[p] = i; + H_csc.x[p] = val; + p++; + } + } + } + H_csc.col_start[n_local] = p; + } + + // Step 2: Factorize H = P * L * D * L^T * P^T + dual_simplex::simplex_solver_settings_t ldlt_settings; + std::vector ldlt_perm; + dual_simplex::csc_matrix_t L_factor(n, n, 1); + std::vector D_factor; + f_t ldlt_work = 0; + f_t ldlt_start = dual_simplex::tic(); + + i_t rank = dual_simplex::right_looking_ldlt( + H_csc, ldlt_settings, f_t(1e-12), ldlt_start, ldlt_perm, L_factor, D_factor, ldlt_work); + + // ldlt_settings uses default time_limit=inf and concurrent_halt=nullptr, + // so only INDEFINITE_MATRIX_RETURN is possible as a negative return code. + cuopt_expects(rank != INDEFINITE_MATRIX_RETURN, error_type_t::ValidationError, - "Quadratic constraint '%s' rotated SOC Q must have at least 1 tail entry", + "Quadratic constraint '%s' is non-convex (Q matrix is indefinite)", qc.constraint_row_name.c_str()); - cone.reserve(q_nnz); - cone.push_back(a); - cone.push_back(b); - cone.insert(cone.end(), tail_vars.begin(), tail_vars.end()); - cone_dim = static_cast(cone.size()); - is_rotated = 1; - rotated_cones.push_back(rotated_soc_t{a, b, tail_vars, false, head_lift_sqrt_ratio}); - } + // Since q_nnz >= 1 is enforced above, Q is nonzero and rank must be >= 1. + // (A nonzero Q entry produces a nonzero H diagonal or off-diagonal, guaranteeing rank > 0.) + assert(rank >= 1); + + // Step 4: Build standard SOC of dimension rank + 2. + // New variables: y_0,...,y_{r-1}, s_0 (head), s_{r+1} (tail) + // Linking rows: + // y_k - sqrt(D[k]) * (L^T P x)_k = 0 for k = 0,...,r-1 + // s_0 + c^T x = alpha + 1/2 + // s_{r+1} + c^T x = alpha - 1/2 + + const i_t r = rank; + const i_t n_new_vars = r + 2; // y_0..y_{r-1}, s_0, s_{r+1} + const i_t n_new_rows = r + 2; + const i_t var_base = csr_A.n; // first new variable index + const i_t y_base = var_base; + const i_t s0_idx = var_base + r; + const i_t sr1_idx = var_base + r + 1; + + // Extend problem dimensions + const f_t pos_inf = std::numeric_limits::infinity(); + const f_t neg_inf = -pos_inf; + user_problem.objective.resize(var_base + n_new_vars, 0); + user_problem.lower.resize(var_base + n_new_vars, neg_inf); + user_problem.upper.resize(var_base + n_new_vars, pos_inf); + user_problem.var_types.resize(var_base + n_new_vars, + dual_simplex::variable_type_t::CONTINUOUS); + if (!user_problem.col_names.empty()) { + user_problem.col_names.resize(var_base + n_new_vars); + for (i_t k = 0; k < r; k++) { + user_problem.col_names[y_base + k] = + "_CUOPT_qc_y_" + std::to_string(qc_i) + "_" + std::to_string(k); + } + user_problem.col_names[s0_idx] = "_CUOPT_qc_s0_" + std::to_string(qc_i); + user_problem.col_names[sr1_idx] = "_CUOPT_qc_sr1_" + std::to_string(qc_i); + } + // s_0 (cone head) — do NOT set lower=0 here; cone membership implies s_0 >= 0 + // and the barrier solver's bound-split logic handles this automatically. - cuopt_expects(have_uniform_s && uniform_s > tol, - error_type_t::ValidationError, - "Quadratic constraint '%s' SOC Q: uniform scale s must be positive (got %.17g)", - qc.constraint_row_name.c_str(), - static_cast(uniform_s)); - qc_soc_uniform_scale[qc_i] = uniform_s; + csr_A.n = var_base + n_new_vars; + is_cone_var.resize(var_base + n_new_vars, 0); - for (const i_t var : cone) { - cuopt_expects(var >= 0 && var < static_cast(is_cone_var.size()), - error_type_t::ValidationError, - "SOC variable index %d is outside [0, %zu)", - static_cast(var), - is_cone_var.size()); + // Extend row storage + const i_t m_before = csr_A.m; + user_problem.rhs.resize(m_before + n_new_rows); + user_problem.row_sense.resize(m_before + n_new_rows); + if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_before + n_new_rows); } + + dual_simplex::sparse_vector_t eq_row; + eq_row.n = csr_A.n; + + // y-linking rows: y_k - sqrt(D[k]) * [row k of L^T P] * x = 0 + // L is unit lower triangular in permuted local indices. + // Column k of L has: L(k,k)=1 at local perm[k], L(j,k) at local perm[j] for j>k. + // (L^T P x)_k = sum_j L(j,k) * x_{var_set[perm[j]]} + for (i_t k = 0; k < r; k++) { + const f_t sqrt_dk = std::sqrt(D_factor[k]); + eq_row.i.clear(); + eq_row.x.clear(); + // y_k coefficient + eq_row.i.push_back(y_base + k); + eq_row.x.push_back(f_t(1)); + // -sqrt(D[k]) * L(:,k) entries applied to x_{var_set[perm[j]]} + for (i_t p = L_factor.col_start[k]; p < L_factor.col_start[k + 1]; p++) { + const i_t j = L_factor.i[p]; // permuted local row index + const f_t l_val = L_factor.x[p]; + const i_t global_var = var_set[ldlt_perm[j]]; + eq_row.i.push_back(global_var); + eq_row.x.push_back(-sqrt_dk * l_val); + } + eq_row.sort(); + csr_A.append_row(eq_row); + user_problem.row_sense[m_before + k] = 'E'; + user_problem.rhs[m_before + k] = 0; + if (!user_problem.row_names.empty()) { + user_problem.row_names[m_before + k] = + "_CUOPT_qc_y_link_" + std::to_string(qc_i) + "_" + std::to_string(k); + } + } + + // s_0 linking row: s_0 + c^T x = alpha + 1/2 + { + eq_row.i.clear(); + eq_row.x.clear(); + eq_row.i.push_back(s0_idx); + eq_row.x.push_back(f_t(1)); + for (size_t p = 0; p < qc.linear_values.size(); ++p) { + const f_t v = qc.linear_values[p]; + if (std::abs(v) < tol) continue; + eq_row.i.push_back(qc.linear_indices[p]); + eq_row.x.push_back(v); + } + eq_row.sort(); + csr_A.append_row(eq_row); + user_problem.row_sense[m_before + r] = 'E'; + user_problem.rhs[m_before + r] = alpha + f_t(0.5); + if (!user_problem.row_names.empty()) { + user_problem.row_names[m_before + r] = "_CUOPT_qc_s0_link_" + std::to_string(qc_i); + } + } + + // s_{r+1} linking row: s_{r+1} + c^T x = alpha - 1/2 + { + eq_row.i.clear(); + eq_row.x.clear(); + eq_row.i.push_back(sr1_idx); + eq_row.x.push_back(f_t(1)); + for (size_t p = 0; p < qc.linear_values.size(); ++p) { + const f_t v = qc.linear_values[p]; + if (std::abs(v) < tol) continue; + eq_row.i.push_back(qc.linear_indices[p]); + eq_row.x.push_back(v); + } + eq_row.sort(); + csr_A.append_row(eq_row); + user_problem.row_sense[m_before + r + 1] = 'E'; + user_problem.rhs[m_before + r + 1] = alpha - f_t(0.5); + if (!user_problem.row_names.empty()) { + user_problem.row_names[m_before + r + 1] = "_CUOPT_qc_sr1_link_" + std::to_string(qc_i); + } + } + + // Register the cone: standard SOC, dim = r+2, head = s_0, tails = (y_0,...,y_{r-1}, s_{r+1}) + cone.clear(); + cone.reserve(r + 2); + cone.push_back(s0_idx); + for (i_t k = 0; k < r; k++) { + cone.push_back(y_base + k); + } + cone.push_back(sr1_idx); + cone_dim = r + 2; + is_rotated = 0; + + for (const i_t var : cone) { + is_cone_var[var] = 1; + } + cone_dims.push_back(cone_dim); + cone_vars.push_back(std::move(cone)); + cone_is_rotated.push_back(is_rotated); } - cone_dims.push_back(cone_dim); - cone_vars.push_back(std::move(cone)); - cone_is_rotated.push_back(is_rotated); } + + // Recount affine linear aux variables (some may have been invalidated by the general path) + n_affine_linear_aux = 0; + for (size_t qc_i = 0; qc_i < qcs.size(); ++qc_i) { + if (qc_affine_heads[qc_i] >= 0) { ++n_affine_linear_aux; } + } + // Add affine linear auxiliary variables and linking rows. if (n_affine_linear_aux > 0) { const f_t inf = std::numeric_limits::infinity(); @@ -528,9 +837,9 @@ void convert_quadratic_constraints_to_second_order_cones( user_problem.row_sense.resize(m_aug); if (!user_problem.row_names.empty()) { user_problem.row_names.resize(m_aug); } - csr_A.n = n_aug; + csr_A.n = std::max(csr_A.n, n_aug); dual_simplex::sparse_vector_t eq_row; - eq_row.n = n_aug; + eq_row.n = csr_A.n; for (size_t qc_i = 0; qc_i < qcs.size(); ++qc_i) { const i_t aux_j = qc_affine_heads[qc_i]; @@ -567,7 +876,7 @@ void convert_quadratic_constraints_to_second_order_cones( "Internal error: CSR row count after affine QC linking"); } - i_t n_prob = n_with_affine_aux; + i_t n_prob = csr_A.n; // Convert rotated SOC cones to standard SOC cones. if (!rotated_cones.empty()) { diff --git a/cpp/src/dual_simplex/right_looking_lu.cpp b/cpp/src/dual_simplex/right_looking_lu.cpp index 34b4ba0ac4..00f2fefb62 100644 --- a/cpp/src/dual_simplex/right_looking_lu.cpp +++ b/cpp/src/dual_simplex/right_looking_lu.cpp @@ -1098,6 +1098,690 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, return pivots; } +// ============================================================================= +// Symmetric positive semidefinite LDL^T factorization with Markowitz pivoting. +// Computes P * A * P^T = L * D * L^T where: +// - A is symmetric PSD (lower triangle stored in CSC) +// - P is a symmetric fill-reducing permutation (from Markowitz) +// - L is unit lower triangular (stored in CSC) +// - D is diagonal (stored as a vector) +// ============================================================================= + +namespace { + +// Represents the lower triangle (including diagonal) of a symmetric trailing matrix +// during right-looking LDL^T factorization. +// Stores column representation of the lower triangle and row representation (index only) +// of the lower triangle. Since the matrix is symmetric, "row i" stores columns j <= i +// that have a nonzero in position (i, j). +template +class symmetric_trailing_matrix_t { + public: + // Construct from a symmetric matrix stored as lower triangle in CSC format. + // A.col_start[j]..A.col_start[j+1]-1 contain entries (i, j) with i >= j. + symmetric_trailing_matrix_t(const csc_matrix_t& A) + : n_(A.n), + Bnz_(0), + work_estimate_(0), + col_start_(n_), + col_end_(n_), + col_max_(n_), + row_start_(n_), + row_end_(n_), + row_max_(n_), + diag_(n_, 0.0), + pivot_col_val_(n_, 0.0), + pivot_col_mark_(n_, 0), + counts_(compute_degree(A), n_), + unused_col_nz_(0), + unused_row_nz_(0) + { + // Count total nonzeros (lower triangle including diagonal) + for (i_t j = 0; j < n_; j++) { + Bnz_ += A.col_start[j + 1] - A.col_start[j]; + } + work_estimate_ += 2 * n_; + + // Allocate 2x initial size for column and row storage + i_t col_nz = 2 * Bnz_; + i_t row_nz = 2 * Bnz_; + + c_i_.resize(col_nz); + c_x_.resize(col_nz); + r_j_.resize(row_nz); + + // Initialize row storage pointers + // Row i stores columns j < i (off-diagonal entries in row i of lower triangle). + // The diagonal is stored separately in diag_. + // Row degree = number of off-diagonal entries in row i of the lower triangle. + i_t nz = 0; + for (i_t i = 0; i < n_; i++) { + row_start_[i] = nz; + row_end_[i] = nz; + // The degree from compute_degree counts off-diag entries in column i (below diagonal) + // plus off-diag entries in row i (same thing by symmetry). + // For row storage, we store columns j < i that appear in column j at row i. + // We'll compute this during column init below. For now, reserve based on degree. + i_t row_space = 2 * counts_.get_count(i); + row_max_[i] = nz + row_space; + nz += row_space; + } + // Resize row storage if needed + if (nz > row_nz) { r_j_.resize(nz); } + work_estimate_ += 4 * n_; + + // Initialize column storage and populate row indices + nz = 0; + for (i_t j = 0; j < n_; j++) { + const i_t A_start = A.col_start[j]; + const i_t A_end = A.col_start[j + 1]; + // Count off-diagonal entries (entries below diagonal) + i_t off_diag_count = 0; + for (i_t p = A_start; p < A_end; p++) { + if (A.i[p] == j) { + diag_[j] = A.x[p]; + } else { + off_diag_count++; + } + } + i_t col_space = 2 * std::max(off_diag_count, i_t(1)); + col_max_[j] = nz + col_space; + col_start_[j] = nz; + col_end_[j] = nz; + + // Store only off-diagonal entries in the column (i > j) + for (i_t p = A_start; p < A_end; p++) { + const i_t row = A.i[p]; + const f_t val = A.x[p]; + if (row == j) { continue; } // diagonal handled above + assert(row > j); // lower triangle: row > col for off-diag + c_i_[col_end_[j]] = row; + c_x_[col_end_[j]] = val; + col_end_[j]++; + + // Also add to row storage: row `row` has an entry in column j + ensure_row_space(row, 1); + r_j_[row_end_[row]] = j; + row_end_[row]++; + } + nz += col_space; + } + // Resize column storage to actual allocated size + if (nz > col_nz) { + c_i_.resize(nz); + c_x_.resize(nz); + } + work_estimate_ += 7 * n_ + 7 * Bnz_; + } + + f_t record_and_clear_work_estimate_() + { + const f_t counts_work = counts_.record_and_clear_work_estimate_(); + work_estimate_ += counts_work; + f_t tmp = work_estimate_; + work_estimate_ = 0; + return tmp; + } + + // Symmetric Markowitz search: find pivot p on diagonal with |diag_[p]| >= pivot_tol + // minimizing degree. Within a fixed degree, select the largest |diag[p]|. + // Degree here is the number of off-diagonal neighbors. + // Returns the number of candidates searched. + i_t symmetric_markowitz_search(f_t pivot_tol, i_t& pivot_p, f_t& pivot_val) + { + i_t best_degree = n_; + i_t nsearch = 0; + + for (i_t nz = 0; nz <= best_degree; nz++) { + const auto& elements = counts_.get_elements_with_count(nz); + for (const i_t p : elements) { + const f_t d = diag_[p]; + if (std::abs(d) >= pivot_tol) { + if (nz < best_degree || (nz == best_degree && std::abs(d) > std::abs(pivot_val))) { + best_degree = nz; + pivot_p = p; + pivot_val = d; + } + } + nsearch++; + } + work_estimate_ += 4 * static_cast(elements.size()); + // Once we've found a pivot at degree nz, no need to search higher degrees + if (pivot_p != -1 && nz == best_degree) { break; } + } + return nsearch; + } + + // Symmetric Schur complement: update A <- A - d * l * l^T (lower triangle only) + // where l_i = A(i, pivot_p) / d for all i != pivot_p, d = diag_[pivot_p]. + // The full pivot vector l comes from: + // - Column pivot_p entries (rows i > pivot_p): stored directly + // - Row pivot_p entries (cols j < pivot_p): these represent A(pivot_p, j) = A(j, pivot_p) + // by symmetry, but are stored in column j at row pivot_p. + void symmetric_schur_complement(i_t pivot_p, f_t drop_tol, f_t pivot_val) + { + // Step 1: Build the full pivot vector from both column (below) and row (left). + // For i > pivot_p: l_i = A(i, pivot_p) / pivot_val, found in column pivot_p + // For j < pivot_p: l_j = A(pivot_p, j) / pivot_val = A(j, pivot_p) / pivot_val + // A(pivot_p, j) is stored in column j at row pivot_p (lower triangle: pivot_p > j) + i_t pivot_col_count = 0; + + // From column storage: entries (i, pivot_p) with i > pivot_p + const i_t c_pivot_start = col_start_[pivot_p]; + const i_t c_pivot_end = col_end_[pivot_p]; + for (i_t p = c_pivot_start; p < c_pivot_end; p++) { + const i_t i = c_i_[p]; + const f_t li = c_x_[p] / pivot_val; + pivot_col_val_[i] = li; + pivot_col_mark_[i] = 1; + pivot_col_index_.push_back(i); + pivot_col_count++; + } + + // From row storage: entries (pivot_p, j) with j < pivot_p + // The value A(pivot_p, j) is stored in column j at row pivot_p. + const i_t r_pivot_start = row_start_[pivot_p]; + const i_t r_pivot_end = row_end_[pivot_p]; + for (i_t rp = r_pivot_start; rp < r_pivot_end; rp++) { + const i_t j = r_j_[rp]; + // Look up A(pivot_p, j) from column j + f_t val = 0; + for (i_t q = col_start_[j]; q < col_end_[j]; q++) { + if (c_i_[q] == pivot_p) { + val = c_x_[q]; + break; + } + } + const f_t lj = val / pivot_val; + pivot_col_val_[j] = lj; + pivot_col_mark_[j] = 1; + pivot_col_index_.push_back(j); + pivot_col_count++; + } + work_estimate_ += 5 * (c_pivot_end - c_pivot_start) + 5 * (r_pivot_end - r_pivot_start); + + // Step 2: For each node j in the pivot vector, update the trailing matrix. + // The update is: A(i, j) -= pivot_val * l_i * l_j for all pairs (i, j) with i >= j + // that are both in the pivot vector. Also update diagonals. + for (i_t k = 0; k < pivot_col_count; k++) { + const i_t j = pivot_col_index_[k]; + const f_t lj = pivot_col_val_[j]; + + // Update diagonal: A(j,j) -= pivot_val * lj * lj + diag_[j] -= pivot_val * lj * lj; + + // Update off-diagonal entries in column j: A(i, j) -= pivot_val * l_i * l_j + // for i > j where l_i != 0 (i.e., i is in pivot vector) + // Column j stores entries with row > j in lower triangle. + + // Count how many pivot vector entries with index > j exist (potential updates/fills) + i_t n_pivot_entries_below_j = 0; + for (i_t m = 0; m < pivot_col_count; m++) { + if (pivot_col_index_[m] > j) { n_pivot_entries_below_j++; } + } + + // Scan existing entries in column j and update those that are in pivot vector + const i_t c_start = col_start_[j]; + const i_t c_end = col_end_[j]; + i_t n_updated = 0; + i_t n_cancel = 0; + for (i_t q = c_start; q < c_end; q++) { + const i_t i = c_i_[q]; + if (pivot_col_mark_[i]) { + // Entry (i, j) exists and i is in pivot vector: update + const f_t li = pivot_col_val_[i]; + c_x_[q] -= pivot_val * li * lj; + if (std::abs(c_x_[q]) < drop_tol) { + c_x_[q] = 0; + n_cancel++; + } + n_updated++; + } + } + work_estimate_ += 4 * (c_end - c_start); + i_t n_fillin = n_pivot_entries_below_j - n_updated; + + // Step 2b: Remove cancellations + if (n_cancel > 0) { + i_t new_end = col_start_[j]; + for (i_t q = col_start_[j]; q < col_end_[j]; q++) { + if (c_x_[q] != 0) { + c_i_[new_end] = c_i_[q]; + c_x_[new_end] = c_x_[q]; + new_end++; + } else { + const i_t dead_row = c_i_[q]; + // Remove column j from row dead_row + for (i_t rp2 = row_start_[dead_row]; rp2 < row_end_[dead_row]; rp2++) { + if (r_j_[rp2] == j) { + r_j_[rp2] = r_j_[row_end_[dead_row] - 1]; + row_end_[dead_row]--; + break; + } + } + // Update degree for dead_row + const i_t rdeg = counts_.get_count(dead_row); + if (rdeg > 0) { counts_.update_count(dead_row, rdeg - 1); } + work_estimate_ += 6; + } + } + col_end_[j] = new_end; + } + + // Step 2c: Insert fill-in entries + if (n_fillin > 0) { + ensure_col_space(j, n_fillin); + for (i_t m = 0; m < pivot_col_count; m++) { + const i_t i = pivot_col_index_[m]; + if (i <= j) { continue; } // only lower triangle: i > j + // Check if (i, j) already exists + bool found = false; + for (i_t q = col_start_[j]; q < col_end_[j]; q++) { + if (c_i_[q] == i) { + found = true; + break; + } + } + if (found) { continue; } + + // Insert fill-in: A(i, j) = -pivot_val * l_i * l_j + const f_t li = pivot_col_val_[i]; + const f_t val = -pivot_val * li * lj; + if (std::abs(val) < drop_tol) { continue; } + + c_i_[col_end_[j]] = i; + c_x_[col_end_[j]] = val; + col_end_[j]++; + + // Insert into row storage: row i gets column j + ensure_row_space(i, 1); + r_j_[row_end_[i]] = j; + row_end_[i]++; + + // Update degree for row i + const i_t rdeg = counts_.get_count(i); + counts_.update_count(i, rdeg + 1); + work_estimate_ += 10; + } + } + + // Update degree for node j + i_t col_entries = col_end_[j] - col_start_[j]; + i_t row_entries = row_end_[j] - row_start_[j]; + i_t total_deg = col_entries + row_entries; + if (total_deg != counts_.get_count(j)) { counts_.update_count(j, total_deg); } + } + + // Step 3: Clear pivot vector workspaces + for (i_t k = 0; k < pivot_col_count; k++) { + const i_t i = pivot_col_index_[k]; + pivot_col_val_[i] = 0; + pivot_col_mark_[i] = 0; + } + pivot_col_index_.clear(); + work_estimate_ += 2 * pivot_col_count; + } + + // Extract the full pivot vector as the L column (divided by pivot value). + // The L column for pivot p includes ALL off-diagonal entries of the full symmetric + // column p: entries below (from column storage) AND entries above (from row storage, + // looked up in other columns). Returns entries in original row indices. + void extract_column(i_t pivot_p, f_t pivot_val, csc_matrix_t& L, i_t& Lnz) + { + // Entries below pivot: stored directly in column pivot_p + const i_t c_start = col_start_[pivot_p]; + const i_t c_end = col_end_[pivot_p]; + for (i_t p = c_start; p < c_end; p++) { + const i_t i = c_i_[p]; + const f_t l_val = c_x_[p] / pivot_val; + L.i.push_back(i); + L.x.push_back(l_val); + Lnz++; + } + + // Entries above pivot: row pivot_p stores columns j < pivot_p. + // The value A(pivot_p, j) is stored in column j at row pivot_p. + const i_t r_start = row_start_[pivot_p]; + const i_t r_end = row_end_[pivot_p]; + for (i_t rp = r_start; rp < r_end; rp++) { + const i_t j = r_j_[rp]; + // Look up A(pivot_p, j) from column j + f_t val = 0; + for (i_t q = col_start_[j]; q < col_end_[j]; q++) { + if (c_i_[q] == pivot_p) { + val = c_x_[q]; + break; + } + } + const f_t l_val = val / pivot_val; + L.i.push_back(j); + L.x.push_back(l_val); + Lnz++; + } + work_estimate_ += 4 * (c_end - c_start) + 5 * (r_end - r_start); + } + + // Remove the pivot node from the trailing matrix: remove column pivot_p and row pivot_p. + void remove_pivot(i_t pivot_p) + { + // Remove pivot_p from all rows that reference it via column storage + // Column pivot_p has entries at rows i > pivot_p. For each such row i, + // remove pivot_p from the row-index list of row i. + const i_t c_start = col_start_[pivot_p]; + const i_t c_end = col_end_[pivot_p]; + for (i_t p = c_start; p < c_end; p++) { + const i_t i = c_i_[p]; + // Remove column pivot_p from row i + for (i_t rp = row_start_[i]; rp < row_end_[i]; rp++) { + if (r_j_[rp] == pivot_p) { + r_j_[rp] = r_j_[row_end_[i] - 1]; + row_end_[i]--; + break; + } + } + // Update degree for row i: decrement by 1 (lost column pivot_p) + const i_t deg = counts_.get_count(i); + if (deg > 0) { counts_.update_count(i, deg - 1); } + } + work_estimate_ += 6 * (c_end - c_start); + + // Remove pivot_p from all columns that row pivot_p references. + // Row pivot_p stores columns j < pivot_p that have an entry at row pivot_p. + const i_t r_start = row_start_[pivot_p]; + const i_t r_end = row_end_[pivot_p]; + for (i_t rp = r_start; rp < r_end; rp++) { + const i_t j = r_j_[rp]; + // Remove row pivot_p from column j + for (i_t q = col_start_[j]; q < col_end_[j]; q++) { + if (c_i_[q] == pivot_p) { + c_i_[q] = c_i_[col_end_[j] - 1]; + c_x_[q] = c_x_[col_end_[j] - 1]; + col_end_[j]--; + break; + } + } + // Update degree for column j: lost row pivot_p + i_t col_entries = col_end_[j] - col_start_[j]; + i_t row_entries = row_end_[j] - row_start_[j]; + i_t total_deg = col_entries + row_entries; + if (total_deg != counts_.get_count(j)) { counts_.update_count(j, total_deg); } + } + work_estimate_ += 6 * (r_end - r_start); + + // Mark pivot as eliminated + col_end_[pivot_p] = col_start_[pivot_p]; + row_end_[pivot_p] = row_start_[pivot_p]; + counts_.remove_from_count(pivot_p); + } + + void garbage_collect(f_t max_unused_fraction = 0.90) + { + if (unused_col_nz_ > max_unused_fraction * static_cast(c_i_.size())) { + std::vector new_c_i; + std::vector new_c_x; + new_c_i.reserve(c_i_.size() - unused_col_nz_); + new_c_x.reserve(c_x_.size() - unused_col_nz_); + for (i_t j = 0; j < n_; j++) { + const i_t new_start = static_cast(new_c_i.size()); + const i_t c_start = col_start_[j]; + const i_t c_end = col_end_[j]; + const i_t col_size = c_end - c_start; + for (i_t p = c_start; p < c_end; p++) { + new_c_i.push_back(c_i_[p]); + new_c_x.push_back(c_x_[p]); + } + col_start_[j] = new_start; + col_end_[j] = static_cast(new_c_i.size()); + for (i_t s = 0; s < col_size; s++) { + new_c_i.push_back(kNone); + new_c_x.push_back(0.0); + } + col_max_[j] = static_cast(new_c_i.size()); + work_estimate_ += 4 * col_size; + } + work_estimate_ += 6 * n_; + c_i_ = std::move(new_c_i); + c_x_ = std::move(new_c_x); + unused_col_nz_ = 0; + } + + if (unused_row_nz_ > max_unused_fraction * static_cast(r_j_.size())) { + std::vector new_r_j; + new_r_j.reserve(r_j_.size() - unused_row_nz_); + for (i_t i = 0; i < n_; i++) { + const i_t new_start = static_cast(new_r_j.size()); + const i_t r_start = row_start_[i]; + const i_t r_end = row_end_[i]; + const i_t row_size = r_end - r_start; + for (i_t p = r_start; p < r_end; p++) { + new_r_j.push_back(r_j_[p]); + } + row_start_[i] = new_start; + row_end_[i] = static_cast(new_r_j.size()); + for (i_t s = 0; s < row_size; s++) { + new_r_j.push_back(kNone); + } + row_max_[i] = static_cast(new_r_j.size()); + work_estimate_ += 2 * row_size; + } + work_estimate_ += 6 * n_; + r_j_ = std::move(new_r_j); + unused_row_nz_ = 0; + } + } + + private: + bool ensure_col_space(i_t j, i_t needed) + { + if (col_end_[j] + needed <= col_max_[j]) { return false; } + const i_t c_start = col_start_[j]; + const i_t c_end = col_end_[j]; + i_t current_size = c_end - c_start; + unused_col_nz_ += current_size; + i_t new_start = c_i_.size(); + for (i_t p = c_start; p < c_end; p++) { + c_i_.push_back(c_i_[p]); + c_x_.push_back(c_x_[p]); + } + work_estimate_ += 2 * (c_end - c_start); + col_start_[j] = new_start; + col_end_[j] = c_i_.size(); + i_t extra = std::max(current_size, needed); + for (i_t k = 0; k < extra; k++) { + c_i_.push_back(kNone); + c_x_.push_back(0.0); + } + work_estimate_ += 2 * extra; + col_max_[j] = c_i_.size(); + work_estimate_ += 10; + return true; + } + + void ensure_row_space(i_t i, i_t needed) + { + if (row_end_[i] + needed <= row_max_[i]) { return; } + const i_t r_start = row_start_[i]; + const i_t r_end = row_end_[i]; + i_t current_size = r_end - r_start; + unused_row_nz_ += current_size; + i_t new_start = r_j_.size(); + for (i_t p = r_start; p < r_end; p++) { + r_j_.push_back(r_j_[p]); + } + work_estimate_ += (r_end - r_start); + row_start_[i] = new_start; + row_end_[i] = r_j_.size(); + i_t extra = std::max(current_size, needed); + for (i_t k = 0; k < extra; k++) { + r_j_.push_back(kNone); + } + work_estimate_ += extra; + row_max_[i] = r_j_.size(); + work_estimate_ += 9; + } + + // Compute the symmetric degree of each node: number of off-diagonal nonzeros + // in the lower triangle that touch node j (entries in column j with row > j, + // plus entries in columns k < j at row j). + std::vector compute_degree(const csc_matrix_t& A) + { + std::vector degree(A.n, 0); + for (i_t j = 0; j < A.n; j++) { + for (i_t p = A.col_start[j]; p < A.col_start[j + 1]; p++) { + const i_t i = A.i[p]; + if (i == j) { continue; } // skip diagonal + // Off-diagonal entry (i, j) with i > j: contributes to degree of both i and j + degree[j]++; + degree[i]++; + } + } + work_estimate_ += 3 * A.n + 2 * Bnz_; + return degree; + } + + i_t n_; + i_t Bnz_; + f_t work_estimate_; + + // Column representation of lower triangle (off-diagonal: rows > col) + std::vector col_start_; + std::vector col_end_; + std::vector col_max_; + std::vector c_i_; + std::vector c_x_; + + // Row representation (index only): row i stores columns j < i with entry (i, j) + std::vector row_start_; + std::vector row_end_; + std::vector row_max_; + std::vector r_j_; + + // Diagonal entries + std::vector diag_; + + // Dense workspaces for pivot column + std::vector pivot_col_val_; + std::vector pivot_col_mark_; + std::vector pivot_col_index_; + + // Single degree structure (symmetric: row degree == col degree) + nonzero_counts_t counts_; + + i_t unused_col_nz_; + i_t unused_row_nz_; +}; + +} // namespace + +template +i_t right_looking_ldlt(const csc_matrix_t& A, + const simplex_solver_settings_t& settings, + f_t pivot_tol, + f_t start_time, + std::vector& perm, + csc_matrix_t& L, + std::vector& D, + f_t& work_estimate) +{ + raft::common::nvtx::range scope("LU::right_looking_ldlt"); + const i_t n = A.n; + assert(A.m == n); + + symmetric_trailing_matrix_t trailing_matrix(A); + + perm.resize(n); + std::fill(perm.begin(), perm.end(), -1); + D.clear(); + D.reserve(n); + + L.m = n; + L.n = n; + L.col_start.resize(n + 1); + L.i.clear(); + L.x.clear(); + + // perminv[original_index] = elimination_order + std::vector perminv(n, -1); + + i_t Lnz = 0; + constexpr f_t drop_tol = 1e-14; + + work_estimate += trailing_matrix.record_and_clear_work_estimate_(); + + i_t pivots = 0; + for (i_t k = 0; k < n; ++k) { + if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) { + return CONCURRENT_HALT_RETURN; + } + if (toc(start_time) > settings.time_limit) { return TIME_LIMIT_RETURN; } + + // Find symmetric pivot + i_t pivot_p = -1; + f_t pivot_val = 0; + trailing_matrix.symmetric_markowitz_search(pivot_tol, pivot_p, pivot_val); + + if (pivot_p == -1) { break; } // No acceptable pivot found (remaining diagonals are zero/tiny) + + // Check for indefiniteness: a negative pivot means the matrix is not PSD + if (pivot_val < 0) { return INDEFINITE_MATRIX_RETURN; } + + // Record permutation + perm[k] = pivot_p; + perminv[pivot_p] = k; + D.push_back(pivot_val); + pivots++; + + // L(:, k) = entries in column pivot_p / pivot_val + L.col_start[k] = Lnz; + // Unit diagonal: L(pivot_p, k) = 1 (implicit, stored for triangular solve compatibility) + L.i.push_back(pivot_p); + L.x.push_back(1.0); + Lnz++; + trailing_matrix.extract_column(pivot_p, pivot_val, L, Lnz); + + // Symmetric Schur complement: A <- A - d * l * l^T + trailing_matrix.symmetric_schur_complement(pivot_p, drop_tol, pivot_val); + + // Remove pivot from trailing matrix + trailing_matrix.remove_pivot(pivot_p); + + trailing_matrix.garbage_collect(); + + work_estimate += trailing_matrix.record_and_clear_work_estimate_(); + } + + // Finalize L + L.col_start[pivots] = Lnz; + // Fill remaining col_start entries for any unpivoted columns + for (i_t k = pivots + 1; k <= n; k++) { + L.col_start[k] = Lnz; + } + + // Complete the permutation for unpivoted nodes + { + i_t next = pivots; + for (i_t i = 0; i < n; i++) { + if (perminv[i] == -1) { + perminv[i] = next; + if (next < n) { perm[next] = i; } + next++; + } + } + } + + // Remap L row indices from original to permuted order + for (i_t p = 0; p < Lnz; p++) { + L.i[p] = perminv[L.i[p]]; + } + work_estimate += 2 * Lnz; + + // Resize L to actual rank + L.n = pivots; + L.nz_max = Lnz; + + return pivots; +} + #ifdef DUAL_SIMPLEX_INSTANTIATE_DOUBLE template int right_looking_lu(const csc_matrix_t& A, @@ -1118,6 +1802,15 @@ template int right_looking_lu_row_permutation_only( double start_time, std::vector& q, std::vector& pinv); + +template int right_looking_ldlt(const csc_matrix_t& A, + const simplex_solver_settings_t& settings, + double pivot_tol, + double start_time, + std::vector& perm, + csc_matrix_t& L, + std::vector& D, + double& work_estimate); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/right_looking_lu.hpp b/cpp/src/dual_simplex/right_looking_lu.hpp index 5f0bf570b8..9b2814dfde 100644 --- a/cpp/src/dual_simplex/right_looking_lu.hpp +++ b/cpp/src/dual_simplex/right_looking_lu.hpp @@ -34,4 +34,26 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, std::vector& q, std::vector& pinv); +// Sparse LDL^T factorization with symmetric Markowitz pivoting. +// Computes P * A * P^T = L * D * L^T for a symmetric positive semidefinite matrix A. +// Input: A is n x n symmetric PSD, stored as lower triangle only in CSC format. +// Output: +// perm[k] = original index of the k-th pivot (length = rank) +// L = unit lower triangular factor (CSC, with permuted row indices) +// D = diagonal factor (length = rank) +// Returns: +// rank >= 0: number of successful pivots with D(k,k) >= pivot_tol (PSD case). +// INDEFINITE_MATRIX_RETURN (-4): a negative pivot was encountered (matrix is not PSD). +// CONCURRENT_HALT_RETURN (-2): concurrent halt requested. +// TIME_LIMIT_RETURN (-3): time limit exceeded. +template +i_t right_looking_ldlt(const csc_matrix_t& A, + const simplex_solver_settings_t& settings, + f_t pivot_tol, + f_t start_time, + std::vector& perm, + csc_matrix_t& L, + std::vector& D, + f_t& work_estimate); + } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/types.hpp b/cpp/src/dual_simplex/types.hpp index ea46a1f67e..b4ff4b8cb3 100644 --- a/cpp/src/dual_simplex/types.hpp +++ b/cpp/src/dual_simplex/types.hpp @@ -23,5 +23,7 @@ constexpr float64_t inf = std::numeric_limits::infinity(); #define CONCURRENT_HALT_RETURN -2 // We return this constant to signal that a time limit has occurred #define TIME_LIMIT_RETURN -3 +// We return this constant to signal that a matrix is indefinite (has a negative pivot) +#define INDEFINITE_MATRIX_RETURN -4 } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/pdlp/solve.cu b/cpp/src/pdlp/solve.cu index 07a4676120..5d29e97fa0 100644 --- a/cpp/src/pdlp/solve.cu +++ b/cpp/src/pdlp/solve.cu @@ -1842,6 +1842,21 @@ optimization_problem_solution_t solve_qcqp( std::get<3>(sol_dual_simplex), std::get<4>(sol_dual_simplex), method_t::Barrier); + + if (has_qc) { + CUOPT_LOG_INFO("Dual variables for problems with quadratic constraints not returned."); + const f_t nan_val = std::numeric_limits::quiet_NaN(); + auto stream = op_problem.get_handle_ptr()->get_stream(); + thrust::fill(rmm::exec_policy(stream), + solution.get_dual_solution().begin(), + solution.get_dual_solution().end(), + nan_val); + thrust::fill(rmm::exec_policy(stream), + solution.get_reduced_cost().begin(), + solution.get_reduced_cost().end(), + nan_val); + } + if (settings.sol_file != "") { CUOPT_LOG_INFO("Writing solution to file %s", settings.sol_file.c_str()); solution.write_to_sol_file(settings.sol_file, op_problem.get_handle_ptr()->get_stream()); diff --git a/cpp/tests/dual_simplex/CMakeLists.txt b/cpp/tests/dual_simplex/CMakeLists.txt index dc4ab35b73..1f6740378f 100644 --- a/cpp/tests/dual_simplex/CMakeLists.txt +++ b/cpp/tests/dual_simplex/CMakeLists.txt @@ -6,4 +6,5 @@ ConfigureTest(DUAL_SIMPLEX_TEST ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/solve.cpp ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/solve_barrier.cu + ${CMAKE_CURRENT_SOURCE_DIR}/unit_tests/right_looking_ldlt.cpp LABELS numopt) diff --git a/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp b/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp new file mode 100644 index 0000000000..7c91a1f72f --- /dev/null +++ b/cpp/tests/dual_simplex/unit_tests/right_looking_ldlt.cpp @@ -0,0 +1,605 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace cuopt::linear_programming::dual_simplex::test { + +// Helper: build a CSC lower-triangle matrix from dense symmetric input (column-major, full matrix). +// Only stores entries (i, j) with i >= j. +static csc_matrix_t dense_to_lower_csc(int n, const std::vector& dense) +{ + // Count nonzeros in lower triangle + int nnz = 0; + for (int j = 0; j < n; j++) { + for (int i = j; i < n; i++) { + if (dense[j * n + i] != 0.0) { nnz++; } + } + } + + csc_matrix_t A(n, n, nnz); + int p = 0; + for (int j = 0; j < n; j++) { + A.col_start[j] = p; + for (int i = j; i < n; i++) { + double val = dense[j * n + i]; + if (val != 0.0) { + A.i[p] = i; + A.x[p] = val; + p++; + } + } + } + A.col_start[n] = p; + return A; +} + +// Helper: verify the factorization P * A * P^T = L * D * L^T +// by computing L * D * L^T and comparing against P * A * P^T. +static void verify_ldlt(int n, + const std::vector& dense_A, + int rank, + const std::vector& perm, + const csc_matrix_t& L, + const std::vector& D, + double tol = 1e-10) +{ + // Build P * A * P^T as dense matrix (rank x rank leading submatrix) + // perm[k] = original index of the k-th pivot + // So P*A*P^T[k, l] = A[perm[k], perm[l]] + std::vector PAPT(rank * rank, 0.0); + for (int k = 0; k < rank; k++) { + for (int l = 0; l < rank; l++) { + PAPT[l * rank + k] = dense_A[perm[l] * n + perm[k]]; + } + } + + // Build L as dense matrix (rank x rank) + // L is stored in CSC with permuted row indices: L.i[p] is the permuted row index. + std::vector L_dense(rank * rank, 0.0); + for (int j = 0; j < rank; j++) { + for (int p = L.col_start[j]; p < L.col_start[j + 1]; p++) { + int i = L.i[p]; + double val = L.x[p]; + if (i < rank) { L_dense[j * rank + i] = val; } + } + } + + // Compute L * D * L^T + std::vector LDLT(rank * rank, 0.0); + for (int i = 0; i < rank; i++) { + for (int j = 0; j < rank; j++) { + double sum = 0.0; + for (int k = 0; k < rank; k++) { + sum += L_dense[k * rank + i] * D[k] * L_dense[k * rank + j]; + } + LDLT[j * rank + i] = sum; + } + } + + // Compare PAPT and LDLT + for (int i = 0; i < rank; i++) { + for (int j = 0; j < rank; j++) { + EXPECT_NEAR(PAPT[j * rank + i], LDLT[j * rank + i], tol) + << "Mismatch at (" << i << ", " << j << ")"; + } + } +} + +// Test 1: 2x2 positive definite diagonal matrix +TEST(right_looking_ldlt, diagonal_2x2) +{ + // A = [4 0; 0 9] + const int n = 2; + std::vector dense = {4.0, 0.0, 0.0, 9.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 2); + EXPECT_EQ(D.size(), 2u); + // Both diagonal entries should appear as pivots + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 2: 3x3 positive definite matrix +TEST(right_looking_ldlt, pd_3x3) +{ + // A = [4 2 1] + // [2 5 3] + // [1 3 6] + const int n = 3; + // Column-major storage of full symmetric matrix + std::vector dense = {4.0, 2.0, 1.0, 2.0, 5.0, 3.0, 1.0, 3.0, 6.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 3); + verify_ldlt(n, dense, rank, perm, L, D); + + // All D entries should be positive (PD matrix) + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } +} + +// Test 3: Rank-1 PSD matrix (singular) +TEST(right_looking_ldlt, rank1_psd) +{ + // A = v * v^T where v = [1, 2, 3] + // A = [1 2 3] + // [2 4 6] + // [3 6 9] + const int n = 3; + std::vector dense = {1.0, 2.0, 3.0, 2.0, 4.0, 6.0, 3.0, 6.0, 9.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + // Rank should be 1 + EXPECT_EQ(rank, 1); + EXPECT_EQ(D.size(), 1u); + EXPECT_GT(D[0], 0.0); + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 4: Rank-2 PSD matrix (singular) +TEST(right_looking_ldlt, rank2_psd) +{ + // A = v1*v1^T + v2*v2^T where v1 = [1, 0, 1], v2 = [0, 1, 1] + // A = [1 0 1] [0 0 0] [1 0 1] + // [0 0 0] + [0 1 1] = [0 1 1] + // [1 0 1] [0 1 1] [1 1 2] + const int n = 3; + std::vector dense = {1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 1.0, 2.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + // Rank should be 2 + EXPECT_EQ(rank, 2); + EXPECT_EQ(D.size(), 2u); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 5: Zero matrix (rank 0) +TEST(right_looking_ldlt, zero_matrix) +{ + const int n = 3; + std::vector dense(n * n, 0.0); + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 0); + EXPECT_TRUE(D.empty()); +} + +// Test 6: 1x1 matrix +TEST(right_looking_ldlt, scalar_1x1) +{ + const int n = 1; + std::vector dense = {7.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 1); + EXPECT_EQ(D.size(), 1u); + EXPECT_NEAR(D[0], 7.0, 1e-14); + EXPECT_EQ(perm[0], 0); +} + +// Test 7: Larger PD matrix (5x5, from A^T*A with random A) +TEST(right_looking_ldlt, pd_5x5) +{ + // Build A = B^T * B where B is 5x5 with known structure + // B = [2 1 0 0 0] + // [1 3 1 0 0] + // [0 1 4 1 0] + // [0 0 1 5 1] + // [0 0 0 1 6] + // A = B^T * B (tridiagonal B -> banded A, guaranteed PD) + const int n = 5; + double B[5][5] = { + {2, 1, 0, 0, 0}, {1, 3, 1, 0, 0}, {0, 1, 4, 1, 0}, {0, 0, 1, 5, 1}, {0, 0, 0, 1, 6}}; + std::vector dense(n * n, 0.0); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + double sum = 0.0; + for (int k = 0; k < n; k++) { + sum += B[k][i] * B[k][j]; + } + dense[j * n + i] = sum; + } + } + + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 5); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 8: Rank-3 PSD 5x5 matrix +TEST(right_looking_ldlt, rank3_5x5_psd) +{ + // A = sum_{k=0}^{2} v_k * v_k^T with linearly independent v_k + // v0 = [1, 0, 1, 0, 1], v1 = [0, 1, 1, 0, 0], v2 = [0, 0, 0, 1, 1] + const int n = 5; + double v[3][5] = {{1, 0, 1, 0, 1}, {0, 1, 1, 0, 0}, {0, 0, 0, 1, 1}}; + std::vector dense(n * n, 0.0); + for (int k = 0; k < 3; k++) { + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + dense[j * n + i] += v[k][i] * v[k][j]; + } + } + } + // A = [1 0 1 0 1] + // [0 1 1 0 0] + // [1 1 2 0 1] + // [0 0 0 1 1] + // [1 0 1 1 2] + + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 3); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 9: Identity matrix +TEST(right_looking_ldlt, identity) +{ + const int n = 4; + std::vector dense(n * n, 0.0); + for (int i = 0; i < n; i++) { + dense[i * n + i] = 1.0; + } + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 4); + for (int k = 0; k < rank; k++) { + EXPECT_NEAR(D[k], 1.0, 1e-14); + } + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 10: Sparse PSD matrix arising from graph Laplacian (always singular, rank = n-1) +TEST(right_looking_ldlt, graph_laplacian) +{ + // 4-node cycle graph Laplacian: + // L = [ 2 -1 0 -1] + // [-1 2 -1 0] + // [ 0 -1 2 -1] + // [-1 0 -1 2] + // Eigenvalues: 0, 2, 2, 4 -> rank 3 + const int n = 4; + std::vector dense = { + 2.0, -1.0, 0.0, -1.0, -1.0, 2.0, -1.0, 0.0, 0.0, -1.0, 2.0, -1.0, -1.0, 0.0, -1.0, 2.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L_out(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L_out, D, work_estimate); + + // Graph Laplacian of connected graph has rank n-1 + EXPECT_EQ(rank, 3); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense, rank, perm, L_out, D); +} + +// Test 11: Verify that unsymmetric input (Q + Q^T) works when symmetrized externally. +// This tests the use case: given arbitrary Q, form H = Q + Q^T (lower triangle) and factorize. +TEST(right_looking_ldlt, symmetrized_from_unsymmetric) +{ + // Q (not symmetric): + // Q = [4 1 0] + // [3 5 2] + // [0 1 6] + // H = Q + Q^T = [8 4 0] + // [4 10 3] + // [0 3 12] + // H is positive definite (diagonally dominant). + const int n = 3; + std::vector dense_H = {8.0, 4.0, 0.0, 4.0, 10.0, 3.0, 0.0, 3.0, 12.0}; + auto A = dense_to_lower_csc(n, dense_H); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, 3); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense_H, rank, perm, L, D); +} + +// Test 12: Larger rank-deficient matrix (10x10, rank 5) +TEST(right_looking_ldlt, rank5_10x10_psd) +{ + const int n = 10; + const int r = 5; + // Build A = V * V^T where V is 10 x 5 + double V[10][5] = {{1, 0, 0, 0, 0}, + {0, 1, 0, 0, 0}, + {0, 0, 1, 0, 0}, + {0, 0, 0, 1, 0}, + {0, 0, 0, 0, 1}, + {1, 1, 0, 0, 0}, + {0, 1, 1, 0, 0}, + {0, 0, 1, 1, 0}, + {0, 0, 0, 1, 1}, + {1, 0, 0, 0, 1}}; + std::vector dense(n * n, 0.0); + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + double sum = 0.0; + for (int k = 0; k < r; k++) { + sum += V[i][k] * V[j][k]; + } + dense[j * n + i] = sum; + } + } + + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank_out = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank_out, r); + for (int k = 0; k < rank_out; k++) { + EXPECT_GT(D[k], 0.0); + } + verify_ldlt(n, dense, rank_out, perm, L, D); +} + +// Test 13: 10x10 matrix that is all zeros except a single 1 at position (9,9) (0-indexed) +TEST(right_looking_ldlt, sparse_single_entry_10x10) +{ + const int n = 10; + std::vector dense(n * n, 0.0); + dense[9 * n + 9] = 1.0; // A(9,9) = 1, everything else is 0 + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + // Only one nonzero diagonal entry, so rank = 1 + EXPECT_EQ(rank, 1); + EXPECT_EQ(D.size(), 1u); + EXPECT_NEAR(D[0], 1.0, 1e-14); + // The pivot should be at original index 9 + EXPECT_EQ(perm[0], 9); + verify_ldlt(n, dense, rank, perm, L, D); +} + +// Test 14: Indefinite matrix (has a negative eigenvalue). +// The factorization should detect this and return INDEFINITE_MATRIX_RETURN. +// Matrix: A = [1, 2; 2, 1] has eigenvalues 3 and -1 (indefinite). +TEST(right_looking_ldlt, indefinite_2x2) +{ + const int n = 2; + std::vector dense = {1.0, 2.0, 2.0, 1.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int result = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + // Should return INDEFINITE_MATRIX_RETURN (-4) since the matrix is indefinite + EXPECT_EQ(result, INDEFINITE_MATRIX_RETURN); +} + +// Test 15: Larger indefinite matrix (4x4 with mixed eigenvalues). +// A = [2, 1, 0, 1; 1, 0, 1, 0; 0, 1, 2, -1; 1, 0, -1, -1] +TEST(right_looking_ldlt, indefinite_4x4) +{ + const int n = 4; + // A symmetric indefinite matrix + std::vector dense = { + 2.0, 1.0, 0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 2.0, -1.0, 1.0, 0.0, -1.0, -1.0}; + auto A = dense_to_lower_csc(n, dense); + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int result = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + // Should return INDEFINITE_MATRIX_RETURN (-4) since the matrix is indefinite + EXPECT_EQ(result, INDEFINITE_MATRIX_RETURN); +} + +// Test 16: Large sparse PD matrix A = e1*e^T + e*e1^T + n*I (n = 10000). +// e1 = [1,0,...,0], e = [1,1,...,1]. +// A(0,0) = n+2, A(i,0) = A(0,i) = 1 for i > 0, A(i,i) = n for i > 0. +// This is an arrowhead matrix: dense first row/column, diagonal elsewhere. +// Markowitz pivoting should eliminate the n-1 diagonal nodes (degree 1) first, +// producing zero fill, then eliminate the dense node last. +// Without reordering, eliminating the dense node first would cause O(n^2) fill. +TEST(right_looking_ldlt, large_arrowhead_markowitz) +{ + const int n = 10000; + + // Build lower triangle CSC of A = e1*e1^T + e*e1^T + n*I + // Column 0: entries at all rows (dense column). A(0,0) = n+2, A(i,0) = 1 for i > 0. + // Columns j > 0: only diagonal entry A(j,j) = n. + // Lower triangle: column 0 has entries (0, n+2), (1, 1), (2, 1), ..., (n-1, 1). + // column j > 0 has entry (j, n). + int nnz = n + (n - 1); // n entries in col 0, (n-1) diagonal entries in cols 1..n-1 + csc_matrix_t A(n, n, nnz); + int p = 0; + // Column 0 + A.col_start[0] = p; + A.i[p] = 0; + A.x[p] = static_cast(n + 2); + p++; + for (int i = 1; i < n; i++) { + A.i[p] = i; + A.x[p] = 1.0; + p++; + } + // Columns 1..n-1 + for (int j = 1; j < n; j++) { + A.col_start[j] = p; + A.i[p] = j; + A.x[p] = static_cast(n); + p++; + } + A.col_start[n] = p; + + simplex_solver_settings_t settings; + std::vector perm; + csc_matrix_t L(n, n, 1); + std::vector D; + double work_estimate = 0; + double start_time = tic(); + + int rank = right_looking_ldlt(A, settings, 1e-12, start_time, perm, L, D, work_estimate); + + EXPECT_EQ(rank, n); + for (int k = 0; k < rank; k++) { + EXPECT_GT(D[k], 0.0); + } + + // With Markowitz pivoting, the diagonal nodes (degree 1) should be eliminated first, + // producing zero fill. The L factor should be very sparse: at most n-1 entries in the + // last column (from the dense node) plus n unit diagonals. + int L_nnz = L.col_start[rank]; + // Each of the first n-1 pivots (degree-1 nodes) produces 1 L entry (the unit diagonal) + // plus 1 off-diagonal entry (connecting to the dense node). The last pivot (dense node) + // produces just the unit diagonal. Total: n diagonals + (n-1) off-diagonals = 2n - 1. + EXPECT_LE(L_nnz, 2 * n); + + // Verify factorization correctness on a few random entries by checking P*A*P^T = L*D*L^T. + // For the arrowhead, after reordering, the first n-1 pivots should be the diagonal nodes. + // Check that perm[0] != 0 (the dense node should NOT be first). + EXPECT_NE(perm[0], 0) << "Markowitz should not pick the dense node (index 0) first"; +} + +} // namespace cuopt::linear_programming::dual_simplex::test diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_test.c b/cpp/tests/linear_programming/c_api_tests/c_api_test.c index 6a587e0da7..d2f73f4073 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_test.c +++ b/cpp/tests/linear_programming/c_api_tests/c_api_test.c @@ -1622,6 +1622,125 @@ cuopt_int_t test_quadratic_constraint_problem(cuopt_int_t* termination_status_pt return status; } +cuopt_int_t test_general_quadratic_constraint_problem(cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr, + cuopt_float_t* solution_values) +{ + cuOptOptimizationProblem problem = NULL; + cuOptSolverSettings settings = NULL; + cuOptSolution solution = NULL; + + // minimize x0 + x1 + // subject to 2*x0^2 + 3*x0*x1 + 2*x1^2 <= 1 (unsymmetric Q, general convex) + // x0 - x1 = 0 + // Q is given with only upper triangle entry for the cross term: + // (0,0,2), (0,1,3), (1,1,2) + // After symmetrization: H = [4 3; 3 4], eigenvalues 1 and 7 (PD). + // With x0 = x1: quadratic = 7*x0^2 <= 1, min 2*x0 at x0 = -1/sqrt(7) + // Optimal objective = -2/sqrt(7) ≈ -0.755929 + cuopt_int_t num_variables = 2; + cuopt_int_t num_linear_constraints = 1; + cuopt_int_t objective_sense = CUOPT_MINIMIZE; + cuopt_float_t objective_offset = 0.0; + cuopt_float_t objective_coefficients[] = {1.0, 1.0}; + + cuopt_int_t row_offsets[] = {0, 2}; + cuopt_int_t column_indices[] = {0, 1}; + cuopt_float_t values[] = {1.0, -1.0}; + + cuopt_float_t constraint_bounds[] = {0.0}; + char constraint_sense[] = {CUOPT_EQUAL}; + + cuopt_float_t var_lower_bounds[] = {-CUOPT_INFINITY, -CUOPT_INFINITY}; + cuopt_float_t var_upper_bounds[] = {CUOPT_INFINITY, CUOPT_INFINITY}; + char variable_types[] = {CUOPT_CONTINUOUS, CUOPT_CONTINUOUS}; + + // Unsymmetric Q: only upper triangle cross term (0,1,3) + cuopt_int_t qc_row_index[] = {0, 0, 1}; + cuopt_int_t qc_col_index[] = {0, 1, 1}; + cuopt_float_t qc_coeff[] = {2.0, 3.0, 2.0}; + + cuopt_int_t status; + + status = cuOptCreateProblem(num_linear_constraints, + num_variables, + objective_sense, + objective_offset, + objective_coefficients, + row_offsets, + column_indices, + values, + constraint_sense, + constraint_bounds, + var_lower_bounds, + var_upper_bounds, + variable_types, + &problem); + + if (status != CUOPT_SUCCESS) { + printf("Error creating problem: %d\n", status); + goto DONE; + } + + status = cuOptAddQuadraticConstraint(problem, + 3, + qc_row_index, + qc_col_index, + qc_coeff, + 0, + NULL, + NULL, + CUOPT_LESS_THAN, + 1.0); + if (status != CUOPT_SUCCESS) { + printf("Error adding quadratic constraint: %d\n", status); + goto DONE; + } + + status = cuOptCreateSolverSettings(&settings); + if (status != CUOPT_SUCCESS) { + printf("Error creating solver settings: %d\n", status); + goto DONE; + } + + status = cuOptSetIntegerParameter(settings, CUOPT_METHOD, CUOPT_METHOD_BARRIER); + if (status != CUOPT_SUCCESS) { + printf("Error setting barrier method: %d\n", status); + goto DONE; + } + + status = cuOptSolve(problem, settings, &solution); + if (status != CUOPT_SUCCESS) { + printf("Error solving problem: %d\n", status); + goto DONE; + } + + status = cuOptGetTerminationStatus(solution, termination_status_ptr); + if (status != CUOPT_SUCCESS) { + printf("Error getting termination status: %d\n", status); + goto DONE; + } + + status = cuOptGetObjectiveValue(solution, objective_ptr); + if (status != CUOPT_SUCCESS) { + printf("Error getting objective value: %d\n", status); + goto DONE; + } + + status = cuOptGetPrimalSolution(solution, solution_values); + if (status != CUOPT_SUCCESS) { + printf("Error getting primal solution: %d\n", status); + goto DONE; + } + +DONE: + cuOptDestroyProblem(&problem); + cuOptDestroySolverSettings(&settings); + cuOptDestroySolution(&solution); + + return status; +} + cuopt_int_t test_write_problem(const char* input_filename, const char* output_filename) { cuOptOptimizationProblem problem = NULL; diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp index 1cf4bbaf9d..9a92f7a310 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp +++ b/cpp/tests/linear_programming/c_api_tests/c_api_tests.cpp @@ -241,6 +241,21 @@ TEST(c_api, test_quadratic_constraint_problem) EXPECT_NEAR(solution_values[3], 5.0, 1e-4); } +TEST(c_api, test_general_quadratic_constraint_problem) +{ + cuopt_int_t termination_status; + cuopt_float_t objective; + cuopt_float_t solution_values[2]; + EXPECT_EQ( + test_general_quadratic_constraint_problem(&termination_status, &objective, solution_values), + CUOPT_SUCCESS); + EXPECT_EQ(termination_status, CUOPT_TERMINATION_STATUS_OPTIMAL); + // Optimal: x0 = x1 = -1/sqrt(7), obj = -2/sqrt(7) ≈ -0.755929 + EXPECT_NEAR(objective, -2.0 / sqrt(7.0), 1e-4); + EXPECT_NEAR(solution_values[0], -1.0 / sqrt(7.0), 1e-4); + EXPECT_NEAR(solution_values[1], -1.0 / sqrt(7.0), 1e-4); +} + TEST(c_api, test_write_problem) { const std::string& rapidsDatasetRootDir = cuopt::test::get_rapids_dataset_root_dir(); diff --git a/cpp/tests/linear_programming/c_api_tests/c_api_tests.h b/cpp/tests/linear_programming/c_api_tests/c_api_tests.h index 4a7ce8dcaf..b80c0ec931 100644 --- a/cpp/tests/linear_programming/c_api_tests/c_api_tests.h +++ b/cpp/tests/linear_programming/c_api_tests/c_api_tests.h @@ -42,6 +42,9 @@ cuopt_int_t test_quadratic_ranged_problem(cuopt_int_t* termination_status_ptr, cuopt_int_t test_quadratic_constraint_problem(cuopt_int_t* termination_status_ptr, cuopt_float_t* objective_ptr, cuopt_float_t* solution_values); +cuopt_int_t test_general_quadratic_constraint_problem(cuopt_int_t* termination_status_ptr, + cuopt_float_t* objective_ptr, + cuopt_float_t* solution_values); cuopt_int_t test_write_problem(const char* input_filename, const char* output_filename); cuopt_int_t test_maximize_problem_dual_variables(cuopt_int_t method, cuopt_int_t* termination_status_ptr, diff --git a/cpp/tests/linear_programming/grpc/grpc_integration_test.cpp b/cpp/tests/linear_programming/grpc/grpc_integration_test.cpp index 59ee0746fc..143b6c76db 100644 --- a/cpp/tests/linear_programming/grpc/grpc_integration_test.cpp +++ b/cpp/tests/linear_programming/grpc/grpc_integration_test.cpp @@ -1242,7 +1242,7 @@ TEST_F(ChunkedUploadTests, ConcurrentChunkedUploads) // unary SubmitJob path) must still round-trip the wire format cleanly and // the worker must surface the SOC validator's ValidationError back to the // client. -TEST_F(ChunkedUploadTests, QuadraticConstraintsUnaryRejectsNonZeroRhs) +TEST_F(ChunkedUploadTests, QuadraticConstraintsUnaryNonZeroRhs) { grpc_client_config_t config; config.timeout_seconds = 60; @@ -1261,15 +1261,10 @@ TEST_F(ChunkedUploadTests, QuadraticConstraintsUnaryRejectsNonZeroRhs) settings.time_limit = 10.0; auto result = client->solve_lp(problem, settings); - // SOC conversion currently requires rhs = 0 on every QC row; QC_Test_1 - // has rhs = 5 / rhs = 10, so the validator rejects it. This proves both - // (a) the QCQP wire format made it intact through the unary submit path - // (otherwise the validator would never have run) and (b) worker error - // propagation correctly forwards the SOC validator's ValidationError to - // the client instead of swallowing it into a fake "successful" response. - EXPECT_FALSE(result.success); - EXPECT_THAT(result.error_message, ::testing::HasSubstr("ValidationError")); - EXPECT_THAT(result.error_message, ::testing::HasSubstr("rhs = 0")); + // QC_Test_1 has rhs = 5 / rhs = 10. The general convex quadratic path + // handles nonzero RHS, so the problem should be accepted and solved. + // This proves the QCQP wire format made it intact through the unary submit path. + EXPECT_TRUE(result.success); } // Force the chunked upload path with both a zero-byte threshold (every array @@ -1289,7 +1284,7 @@ TEST_F(ChunkedUploadTests, QuadraticConstraintsUnaryRejectsNonZeroRhs) // that drops or duplicates QC array bytes would manifest as a *different* // failure mode (typically a malformed-problem error or a successful solve // of a tampered problem) rather than the expected rhs=0 rejection. -TEST_F(ChunkedUploadTests, QuadraticConstraintsChunkedRejectsNonZeroRhs) +TEST_F(ChunkedUploadTests, QuadraticConstraintsChunkedNonZeroRhs) { grpc_client_config_t config; config.timeout_seconds = 60; @@ -1309,9 +1304,9 @@ TEST_F(ChunkedUploadTests, QuadraticConstraintsChunkedRejectsNonZeroRhs) settings.time_limit = 10.0; auto result = client->solve_lp(problem, settings); - EXPECT_FALSE(result.success); - EXPECT_THAT(result.error_message, ::testing::HasSubstr("ValidationError")); - EXPECT_THAT(result.error_message, ::testing::HasSubstr("rhs = 0")); + // QC_Test_1 has nonzero RHS, now handled by the general convex quadratic path. + // This proves the chunked wire format correctly transmits QC data. + EXPECT_TRUE(result.success); } // End-to-end SOCP correctness via gRPC: QC_Test_2 is a small convex QCQP diff --git a/cpp/tests/socp/CMakeLists.txt b/cpp/tests/socp/CMakeLists.txt index d53049b2d2..f380984225 100644 --- a/cpp/tests/socp/CMakeLists.txt +++ b/cpp/tests/socp/CMakeLists.txt @@ -6,4 +6,5 @@ ConfigureTest(SOCP_TEST ${CMAKE_CURRENT_SOURCE_DIR}/second_order_cone_kernels.cu ${CMAKE_CURRENT_SOURCE_DIR}/solve_barrier_socp.cu + ${CMAKE_CURRENT_SOURCE_DIR}/general_quadratic_test.cu LABELS numopt) diff --git a/cpp/tests/socp/general_quadratic_test.cu b/cpp/tests/socp/general_quadratic_test.cu new file mode 100644 index 0000000000..2918d80a7d --- /dev/null +++ b/cpp/tests/socp/general_quadratic_test.cu @@ -0,0 +1,827 @@ +/* clang-format off */ +/* + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION & AFFILIATES. All rights reserved. + * SPDX-License-Identifier: Apache-2.0 + */ +/* clang-format on */ + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +namespace cuopt::linear_programming::detail::test { + +using i_t = int; +using f_t = double; +using qc_t = optimization_problem_interface_t::quadratic_constraint_t; + +static void init_handler(const raft::handle_t* handle_ptr) +{ + RAFT_CUBLAS_TRY(raft::linalg::detail::cublassetpointermode( + handle_ptr->get_cublas_handle(), CUBLAS_POINTER_MODE_DEVICE, handle_ptr->get_stream())); + RAFT_CUSPARSE_TRY(raft::sparse::detail::cusparsesetpointermode( + handle_ptr->get_cusparse_handle(), CUSPARSE_POINTER_MODE_DEVICE, handle_ptr->get_stream())); +} + +// Test: general convex quadratic constraint with dense PD Q matrix. +// minimize x0 + x1 +// subject to x^T Q x <= 1 where Q = [2 1; 1 2] (PD, eigenvalues 1 and 3) +// The feasible set is an ellipse centered at origin. +// Optimal should be at the boundary of the ellipse along direction (-1, -1). +// For Q = [2 1; 1 2], the minimum of (x0+x1) on x^T Q x <= 1 is: +// c = (1,1), Q^{-1} c = (1/3)(1,1), c^T Q^{-1} c = 2/3 +// min c^T x = -sqrt(c^T Q^{-1} c) = -sqrt(2/3) +TEST(general_quadratic, dense_pd_2x2_solve) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + // Need at least one linear constraint for the barrier solver. + // Use equality x0 - x1 = 0 to force x0 = x1. + // With x0=x1=t: x^T[2 1;1 2]x = 6t^2 <= 1, obj = 2t. + // Optimal: t = -1/sqrt(6), obj = -2/sqrt(6) = -sqrt(2/3). + constexpr int m = 1; + constexpr int n = 2; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 1.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = -1.0; + + user_problem.rhs = {0.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Build quadratic constraint: x^T [2 1; 1 2] x <= 1 + // Q in COO (lower triangular stored): + // (0,0,2), (1,0,1), (1,1,2) + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "ellipse"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 1.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {2.0, 1.0, 2.0}; + + // Convert to CSR for translation (must include the linear constraint row) + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, -1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + + // Convert CSR back to CSC for the barrier solver + csr_A.to_compressed_col(user_problem.A); + + // Verify a cone was created + EXPECT_GT(user_problem.second_order_cone_dims.size(), 0u); + EXPECT_GT(user_problem.cone_var_start, 0); + + // Debug: verify user_problem dimensions are consistent + EXPECT_EQ(user_problem.A.m, user_problem.num_rows); + EXPECT_EQ(user_problem.A.n, user_problem.num_cols); + EXPECT_EQ(static_cast(user_problem.objective.size()), user_problem.num_cols); + EXPECT_EQ(static_cast(user_problem.lower.size()), user_problem.num_cols); + EXPECT_EQ(static_cast(user_problem.upper.size()), user_problem.num_cols); + EXPECT_EQ(static_cast(user_problem.var_types.size()), user_problem.num_cols); + EXPECT_EQ(static_cast(user_problem.rhs.size()), user_problem.num_rows); + EXPECT_EQ(static_cast(user_problem.row_sense.size()), user_problem.num_rows); + + // Verify cone layout: cone vars should be a trailing block + i_t cone_end = user_problem.cone_var_start; + for (i_t d : user_problem.second_order_cone_dims) { + cone_end += d; + } + EXPECT_EQ(cone_end, user_problem.num_cols) + << "cone_var_start=" << user_problem.cone_var_start + << " cone_dims sum=" << (cone_end - user_problem.cone_var_start) + << " num_cols=" << user_problem.num_cols; + + // Debug: verify row senses are all 'E' (no inequality rows that would generate slacks) + int n_L_rows = 0; + for (int i = 0; i < user_problem.num_rows; ++i) { + if (user_problem.row_sense[i] == 'L') n_L_rows++; + } + EXPECT_EQ(n_L_rows, 0) << "Expected all rows to be equality after conversion, but found " + << n_L_rows << " 'L' rows out of " << user_problem.num_rows; + + // Now solve via barrier + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + // min (x0+x1) s.t. (1/2)*x^T*[4,1;1,4]*x <= 1 with x0=x1 + // With x0=x1=t: (1/2)*(4t^2+t^2+t^2+4t^2) = 5t^2 <= 1 + // Min 2t at t = -1/sqrt(5), obj = -2/sqrt(5) + EXPECT_NEAR(solution.objective, -2.0 / std::sqrt(5.0), 1e-4); +} + +// Test: non-convex quadratic constraint should be rejected. +// Q = [1 2; 2 1] has eigenvalues 3 and -1 (indefinite). +TEST(general_quadratic, rejects_non_convex) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 0; + constexpr int n = 2; + constexpr int nz = 0; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 0}; + + user_problem.rhs.clear(); + user_problem.row_sense.clear(); + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: (0,0,1), (1,0,4), (1,1,1) → H(0,0)=2, H(1,0)=4, H(1,1)=2 + // Full H = [2 4; 4 2], eigenvalues 6 and -2 → indefinite + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "non_convex"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 1.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {1.0, 4.0, 1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0}; + + std::vector qcs = {qc}; + + // Should throw validation error for non-convex Q + EXPECT_THROW( + (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem)), + cuopt::logic_error); +} + +// Test: rank-deficient PSD Q (e.g., Q = v*v^T with v = [1, 1]) +// minimize x0 + x1 +// subject to (x0 + x1)^2 <= 4 (i.e., |x0 + x1| <= 2) +// Q = [1 1; 1 1] has rank 1, eigenvalues 0 and 2. +// Optimal: x0 + x1 = -2, objective = -2 +TEST(general_quadratic, rank_deficient_psd_solve) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 2; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = -1.0; + + user_problem.rhs = {0.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q = [1 2; 2 1] gives H = [2 2; 2 2] (rank 1), rhs = 4 + // (1/2)*x^T*[2,2;2,2]*x = (x0+x1)^2 <= 4 + // With x0=x1=t: (2t)^2 <= 4 → t >= -1, obj = 2t = -2 + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rank1_cone"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 4.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {1.0, 2.0, 1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, -1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + + // Convert CSR back to CSC for the barrier solver + csr_A.to_compressed_col(user_problem.A); + + // The cone dimension should be rank + 2 = 1 + 2 = 3 + ASSERT_EQ(user_problem.second_order_cone_dims.size(), 1u); + EXPECT_EQ(user_problem.second_order_cone_dims[0], 3); + + // Solve + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + // x0=x1=t from equality. Objective = x0 = t. + // Quadratic form from COO (0,0,1),(1,0,2),(1,1,1): H=[2,2;2,2]. + // (1/2)*x^T*H*x = (1/2)*(2t^2+2t^2+2t^2+2t^2) = 4t^2 <= 4, so t >= -1. + // min x0 = min t = -1. + EXPECT_NEAR(solution.objective, -1.0, 1e-4); +} + +// Test: general quadratic constraint WITH an inequality linear constraint. +// minimize x0 + x1 +// subject to x^T [2 1; 1 2] x <= 1 (quadratic, via general path) +// x0 + x1 <= 10 (linear inequality) +// x0 - x1 = 0 (linear equality) +// The inequality is inactive at optimum, so same answer as dense_pd_2x2_solve. +TEST(general_quadratic, with_inequality_constraint) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + // 2 constraints: one equality, one inequality + constexpr int m = 2; + constexpr int n = 2; + constexpr int nz = 4; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 1.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + // Col 0: rows 0 and 1. Col 1: rows 0 and 1. + user_problem.A.col_start = {0, 2, 4}; + user_problem.A.i[0] = 0; // row 0: x0 - x1 = 0 + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 1; // row 1: x0 + x1 <= 10 + user_problem.A.x[1] = 1.0; + user_problem.A.i[2] = 0; // row 0: x0 - x1 = 0 + user_problem.A.x[2] = -1.0; + user_problem.A.i[3] = 1; // row 1: x0 + x1 <= 10 + user_problem.A.x[3] = 1.0; + + user_problem.rhs = {0.0, 10.0}; + user_problem.row_sense = {'E', 'L'}; + user_problem.lower = {-inf, -inf}; + user_problem.upper = {inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: (0,0,2), (1,0,1), (1,1,2) → H = [4,1;1,4] + // (1/2) x^T [4,1;1,4] x <= 1 + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "ellipse_ineq"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 1.0; + qc.rows = {0, 1, 1}; + qc.cols = {0, 0, 1}; + qc.vals = {2.0, 1.0, 2.0}; + + // Build CSR matching the A matrix + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2, 4}; + csr_A.j = {0, 1, 0, 1}; + csr_A.x = {1.0, -1.0, 1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + + // Convert CSR back to CSC for the barrier solver + csr_A.to_compressed_col(user_problem.A); + + // Verify cone layout + EXPECT_GT(user_problem.second_order_cone_dims.size(), 0u); + i_t cone_end = user_problem.cone_var_start; + for (i_t d : user_problem.second_order_cone_dims) { + cone_end += d; + } + EXPECT_EQ(cone_end, user_problem.num_cols) + << "Cone must be trailing block: cone_var_start=" << user_problem.cone_var_start + << " cone_end=" << cone_end << " num_cols=" << user_problem.num_cols; + + // Solve + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + // Same as dense_pd_2x2 test: -2/sqrt(5) + EXPECT_NEAR(solution.objective, -2.0 / std::sqrt(5.0), 1e-4); +} + +// Test: minimize t subject to ||A*x - b||^2 <= t, with b = A*e (b in range of A). +// Since b is achievable, optimal t* = 0 (at x = e). +// +// A = [1 1; 1 -1; 0 1] (3x2), e = [1, 1], b = A*e = [2, 0, 1] +// A^T*A = [2 0; 0 3], A^T*b = [2, 3], b^T*b = 5 +// Variables: z = (x0, x1, t). Objective: min t → c_obj = (0, 0, 1). +// Quadratic constraint: x^T*(A^T*A)*x - 2*(A^T*b)^T*x - t <= -b^T*b +// Q COO (3x3, only x-block): (0,0,2), (1,1,3) +// linear: (-4, -6, -1) +// rhs: -5 +TEST(general_quadratic, least_squares_b_in_range) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + // Variables: x0, x1, u (u = t - b^T*b = t - 5). + // Linear constraint: x0 + x1 = 2 (one row of Ax = b, helps bound x) + constexpr int m = 1; + constexpr int n = 3; // x0, x1, u + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 0.0, 1.0}; // minimize u (obj = u + 5 = t) + user_problem.obj_constant = 5.0; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {2.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, -5.0}; // u >= -5 (since u = t-5, t >= 0) + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Quadratic constraint: x^T*(A^T*A)*x - 2*(A^T*b)^T*x - u <= 0 + // This is equivalent to ||Ax - b||^2 <= t (with t = s + b^T*b, but we use t directly). + // Reformulated with rhs = 0: x^T*Q*x + c^T*z <= 0 + // Q COO: (0,0,2), (1,1,3) — diagonal of A^T*A + // linear: (-4, -6, -1) — (-2*A^T*b on x, -1 on t) + // rhs: 0 + // Note: ||Ax-b||^2 = x^T*Q*x - 2*(A^T*b)^T*x + b^T*b, so constraint is + // ||Ax-b||^2 - t <= 0 ⟺ x^T*Q*x - 2*(A^T*b)^T*x + b^T*b - t <= 0 + // We absorb b^T*b into the objective: min t, with t = s + b^T*b where s is our variable. + // Simpler: just use variables (x0, x1, s) where s = t - b^T*b, minimize s (obj* = -b^T*b + t*). + // Actually simplest: since b=A*e, optimal is t*=0, just verify obj near 0. + // Let's use rhs = 0 formulation: x^T*Q*x - 2*(A^T*b)^T*x + 5 - t <= 0 + // i.e. x^T*Q*x - 2*(A^T*b)^T*x - t <= -5 ... that's negative rhs again. + // + // Alternative: use variable substitution. Let u = t - 5. Then constraint: + // x^T*Q*x - 4*x0 - 6*x1 - u <= 0 (rhs = 0), and objective = min(u + 5). + // At optimum x=(1,1), constraint: 2+3-4-6-u<=0 → -5-u<=0 → u>=-5, so min u=-5, obj=0. + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "least_squares"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1}; + qc.cols = {0, 1}; + qc.vals = {2.0, 3.0}; + qc.linear_values = {-4.0, -6.0, -1.0}; + qc.linear_indices = {0, 1, 2}; + + // Build CSR with the linear constraint: x0 + x1 = 2 + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + // Verify cone layout + i_t cone_end_ls = user_problem.cone_var_start; + for (i_t d : user_problem.second_order_cone_dims) { + cone_end_ls += d; + } + EXPECT_EQ(cone_end_ls, user_problem.num_cols) + << "cone_var_start=" << user_problem.cone_var_start << " cone_end=" << cone_end_ls + << " num_cols=" << user_problem.num_cols; + + // Check that cone variables have valid bounds for barrier + for (i_t j = user_problem.cone_var_start; j < user_problem.num_cols; ++j) { + EXPECT_TRUE(user_problem.lower[j] == 0.0 || user_problem.lower[j] <= -1e30) + << "cone var " << j << " has invalid lower=" << user_problem.lower[j]; + EXPECT_TRUE(user_problem.upper[j] >= 1e30) + << "cone var " << j << " has invalid upper=" << user_problem.upper[j]; + } + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + // b is in range(A) so optimal t* = 0. + // solution.objective = u* = -5; total objective = u* + obj_constant = -5 + 5 = 0 + EXPECT_NEAR(solution.objective + user_problem.obj_constant, 0.0, 1e-3); +} + +// Test: minimize t subject to ||A*x - b||^2 <= t, with b NOT in range(A). +// Optimal t* = ||A*x* - b||^2 > 0 where x* is the least-squares solution. +// +// A = [1 0; 0 1; 0 0] (3x2), b = [1, 1, 1] +// A^T*A = I_2, A^T*b = [1, 1], b^T*b = 3 +// Least-squares solution: x* = A^T*b = [1, 1], residual = b - A*x* = [0, 0, 1] +// Optimal t* = ||residual||^2 = 1 +// +// Variables: z = (x0, x1, t). Objective: min t. +// Q COO: (0,0,1), (1,1,1) — identity A^T*A +// linear: (-2, -2, -1) — i.e. -2*A^T*b and -1 for t +// rhs: -3 — i.e. -b^T*b +TEST(general_quadratic, least_squares_b_not_in_range) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 3; // x0, x1, t + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {0.0, 0.0, 1.0}; // minimize t + + // Equality: x0 + x1 = 2 (bounds x, but the LS solution x*=[1,1] satisfies this) + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 1, 2, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + user_problem.A.i[1] = 0; + user_problem.A.x[1] = 1.0; + + user_problem.rhs = {2.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, 0.0}; // t >= 0 + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Quadratic constraint: x^T*I*x - 2*[1,1]*x - t <= -3 + // i.e. x0^2 + x1^2 - 2*x0 - 2*x1 - t <= -3 + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "least_squares_nofit"; + qc.constraint_row_type = 'L'; + qc.rhs_value = -3.0; + qc.rows = {0, 1}; + qc.cols = {0, 1}; + qc.vals = {1.0, 1.0}; + qc.linear_values = {-2.0, -2.0, -1.0}; + qc.linear_indices = {0, 1, 2}; + + // Build CSR: x0 + x1 = 2 + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 2}; + csr_A.j = {0, 1}; + csr_A.x = {1.0, 1.0}; + + std::vector qcs = {qc}; + convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem); + csr_A.to_compressed_col(user_problem.A); + + simplex_solver_settings_t settings; + settings.barrier = true; + settings.barrier_presolve = true; + settings.dualize = 0; + + lp_solution_t solution(user_problem.num_rows, user_problem.num_cols); + auto status = solve_linear_program_with_barrier(user_problem, settings, solution); + + EXPECT_EQ(status, lp_status_t::OPTIMAL); + // b is NOT in range(A), optimal t* = ||residual||^2 = 1 + EXPECT_NEAR(solution.objective, 1.0, 1e-3); +} + +// Test: x0^2 + x1^2 - t^2 <= 0 with t >= 0 should be accepted (valid SOC: ||x|| <= t). +TEST(general_quadratic, soc_head_nonneg_accepted) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + // Variables: x0, x1, t. Constraint: x0^2 + x1^2 - t^2 <= 0 + constexpr int m = 1; + constexpr int n = 3; + constexpr int nz = 2; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0, 0.0}; // minimize x0 + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 0, 2}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; // dummy row for barrier: t <= 10 + user_problem.A.i[1] = 0; + user_problem.A.x[1] = 0.0; // placeholder + + // Actually just use: x1 = 1 as a simple equality + user_problem.A.col_start = {0, 0, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + + user_problem.rhs = {1.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, 0.0}; // t >= 0 + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: x0^2 + x1^2 - t^2 <= 0 + // (0,0,1), (1,1,1), (2,2,-1) + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "soc_valid"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2}; + qc.cols = {0, 1, 2}; + qc.vals = {1.0, 1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, 1); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1}; + csr_A.j = {1}; + csr_A.x = {1.0}; + + std::vector qcs = {qc}; + // Should NOT throw — head variable t has lower >= 0 + EXPECT_NO_THROW( + (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem))); + + // Verify it produced a cone + EXPECT_GT(user_problem.second_order_cone_dims.size(), 0u); +} + +// Test: x0^2 + x1^2 - t^2 <= 0 with t free should be rejected (non-convex without t >= 0). +TEST(general_quadratic, soc_head_free_rejected) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 3; + constexpr int nz = 1; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + + user_problem.rhs = {1.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, -inf}; // t is FREE — no lower bound + user_problem.upper = {inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: x0^2 + x1^2 - t^2 <= 0 (same Q, but t is free) + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "soc_invalid"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2}; + qc.cols = {0, 1, 2}; + qc.vals = {1.0, 1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1}; + csr_A.j = {1}; + csr_A.x = {1.0}; + + std::vector qcs = {qc}; + // Head variable t is free with no constraint implying t >= 0, so this is non-convex. + EXPECT_THROW( + (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem)), + cuopt::logic_error); +} + +// Test: x0^2 + x1^2 - 2*y*z <= 0 with y >= 0, z >= 0 should be accepted (valid rotated SOC). +TEST(general_quadratic, rotated_soc_heads_nonneg_accepted) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + // Variables: x0, x1, y, z. Constraint: x0^2 + x1^2 - 2*y*z <= 0 + constexpr int m = 1; + constexpr int n = 4; + constexpr int nz = 1; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0, 0.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 1, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + + user_problem.rhs = {1.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, 0.0, 0.0}; // y >= 0, z >= 0 + user_problem.upper = {inf, inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: x0^2 + x1^2 - 2*y*z <= 0 + // Diagonal: (0,0,1), (1,1,1). Off-diagonal: (2,3,-1), (3,2,-1) + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rsoc_valid"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2, 3}; + qc.cols = {0, 1, 3, 2}; + qc.vals = {1.0, 1.0, -1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1}; + csr_A.j = {1}; + csr_A.x = {1.0}; + + std::vector qcs = {qc}; + // Should NOT throw — both head variables y and z have lower >= 0 + EXPECT_NO_THROW( + (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem))); + + EXPECT_GT(user_problem.second_order_cone_dims.size(), 0u); +} + +// Test: x0^2 + x1^2 - 2*y*z <= 0 with y free, z free should be rejected (non-convex). +TEST(general_quadratic, rotated_soc_heads_free_rejected) +{ + raft::handle_t handle{}; + init_handler(&handle); + + using namespace cuopt::linear_programming::dual_simplex; + user_problem_t user_problem(&handle); + + constexpr int m = 1; + constexpr int n = 4; + constexpr int nz = 1; + + user_problem.num_rows = m; + user_problem.num_cols = n; + user_problem.objective = {1.0, 0.0, 0.0, 0.0}; + + user_problem.A.m = m; + user_problem.A.n = n; + user_problem.A.nz_max = nz; + user_problem.A.reallocate(nz); + user_problem.A.col_start = {0, 0, 1, 1, 1}; + user_problem.A.i[0] = 0; + user_problem.A.x[0] = 1.0; + + user_problem.rhs = {1.0}; + user_problem.row_sense = {'E'}; + user_problem.lower = {-inf, -inf, -inf, -inf}; // y and z are FREE + user_problem.upper = {inf, inf, inf, inf}; + user_problem.num_range_rows = 0; + user_problem.var_types.assign(n, variable_type_t::CONTINUOUS); + + // Q COO: same as above + qc_t qc; + qc.constraint_row_index = 0; + qc.constraint_row_name = "rsoc_invalid"; + qc.constraint_row_type = 'L'; + qc.rhs_value = 0.0; + qc.rows = {0, 1, 2, 3}; + qc.cols = {0, 1, 3, 2}; + qc.vals = {1.0, 1.0, -1.0, -1.0}; + + dual_simplex::csr_matrix_t csr_A(m, n, nz); + csr_A.m = m; + csr_A.n = n; + csr_A.row_start = {0, 1}; + csr_A.j = {1}; + csr_A.x = {1.0}; + + std::vector qcs = {qc}; + // Head variables y and z are free with no constraints implying non-negativity. + EXPECT_THROW( + (convert_quadratic_constraints_to_second_order_cones(n, qcs, csr_A, user_problem)), + cuopt::logic_error); +} + +} // namespace cuopt::linear_programming::detail::test diff --git a/python/cuopt/cuopt/tests/socp/test_socp.py b/python/cuopt/cuopt/tests/socp/test_socp.py index a406f1a2e6..3926114e4c 100644 --- a/python/cuopt/cuopt/tests/socp/test_socp.py +++ b/python/cuopt/cuopt/tests/socp/test_socp.py @@ -165,3 +165,34 @@ def test_socp_3_barrier_solution(): assert h1.Value == pytest.approx(1.0, abs=PRIMAL_TOL) assert h2.Value == pytest.approx(1.0, abs=PRIMAL_TOL) assert h3.Value == pytest.approx(1.0, abs=PRIMAL_TOL) + + +def test_general_quadratic_unsymmetric(): + """ + Min x0 + x1 + s.t. 2*x0^2 + 3*x0*x1 + 2*x1^2 <= 1 (unsymmetric Q: cross term only as x0*x1) + x0 - x1 = 0 + + Q is given unsymmetrically: the 3*x0*x1 term is stored as a single + entry (row=0, col=1, val=3) rather than symmetric (0,1,1.5)+(1,0,1.5). + After symmetrization H = [4 3; 3 4], eigenvalues 1 and 7 (PD). + + With x0 = x1 = t: 2t^2 + 3t^2 + 2t^2 = 7t^2 <= 1 + min 2t at t = -1/sqrt(7), obj = -2/sqrt(7) ≈ -0.755929 + """ + problem = Problem("general_qc_unsymmetric") + x0 = problem.addVariable(lb=-np.inf, name="x0") + x1 = problem.addVariable(lb=-np.inf, name="x1") + problem.setObjective(x0 + x1) + problem.addConstraint(2 * x0 * x0 + 3 * x0 * x1 + 2 * x1 * x1 <= 1) + problem.addConstraint(x0 - x1 == 0) + + solution = _solve(problem) + _assert_solution_on_original_model(problem, solution) + _assert_feasible(problem) + + expected_obj = -2.0 / np.sqrt(7.0) + expected_x = -1.0 / np.sqrt(7.0) + assert problem.ObjValue == pytest.approx(expected_obj, abs=OBJ_TOL) + assert x0.Value == pytest.approx(expected_x, abs=PRIMAL_TOL) + assert x1.Value == pytest.approx(expected_x, abs=PRIMAL_TOL)