diff --git a/.github/julia/runtests.jl b/.github/julia/runtests.jl index dda903ed..6315453c 100644 --- a/.github/julia/runtests.jl +++ b/.github/julia/runtests.jl @@ -18,11 +18,11 @@ import Uno_jll Create a new `AmplNLWriter.Optimizer` object that uses Uno as the backing solver. """ -function Optimizer(options = String["logger=SILENT"]) +function Optimizer(options = String["logger=SILENT", "unbounded_objective_threshold=-1e15"]) return AmplNLWriter.Optimizer(Uno_jll.amplexe, options) end -Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000"]) +Optimizer_LP() = Optimizer(["logger=SILENT", "preset=filterslp", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) # This testset runs https://github.com/jump-dev/MINLPTests.jl @testset "MINLPTests" begin diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index 06a07436..1a58b3c7 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp @@ -3,12 +3,13 @@ #include #include "PrimalDualInteriorPointSubproblem.hpp" -#include "optimization/Direction.hpp" -#include "optimization/Iterate.hpp" #include "ingredients/hessian_models/HessianModelFactory.hpp" #include "linear_algebra/SparseStorageFactory.hpp" +#include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp" #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "solvers/SymmetricIndefiniteLinearSolverFactory.hpp" +#include "optimization/Direction.hpp" +#include "optimization/Iterate.hpp" #include "optimization/WarmstartInformation.hpp" #include "preprocessing/Preprocessing.hpp" #include "reformulation/l1RelaxedProblem.hpp" @@ -123,26 +124,6 @@ namespace uno { void PrimalDualInteriorPointSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { - // barrier Lagrangian Hessian - if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { - // original Lagrangian Hessian - this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints); - - // diagonal barrier terms (grouped by variable) - for (size_t variable_index: Range(problem.number_variables)) { - double diagonal_barrier_term = 0.; - if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded - const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index); - diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound; - } - if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded - const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index); - diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound; - } - this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index); - } - } - // barrier objective gradient if (warmstart_information.objective_changed) { // original objective gradient @@ -174,6 +155,26 @@ namespace uno { problem.evaluate_constraints(current_iterate, this->constraints); problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian); } + + // barrier Lagrangian Hessian + if (warmstart_information.objective_changed || warmstart_information.constraints_changed) { + // original Lagrangian Hessian + this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints); + + // diagonal barrier terms (grouped by variable) + for (size_t variable_index: Range(problem.number_variables)) { + double diagonal_barrier_term = 0.; + if (is_finite(problem.variable_lower_bound(variable_index))) { // lower bounded + const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_lower_bound(variable_index); + diagonal_barrier_term += current_multipliers.lower_bounds[variable_index] / distance_to_bound; + } + if (is_finite(problem.variable_upper_bound(variable_index))) { // upper bounded + const double distance_to_bound = current_iterate.primals[variable_index] - problem.variable_upper_bound(variable_index); + diagonal_barrier_term += current_multipliers.upper_bounds[variable_index] / distance_to_bound; + } + this->hessian_model->hessian.insert(diagonal_barrier_term, variable_index, variable_index); + } + } } void PrimalDualInteriorPointSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, @@ -199,7 +200,7 @@ namespace uno { this->evaluate_functions(statistics, problem, current_iterate, current_multipliers, warmstart_information); // compute the primal-dual solution - this->assemble_augmented_system(statistics, problem, current_multipliers); + this->assemble_augmented_system(statistics, problem, current_multipliers, warmstart_information); this->augmented_system.solve(*this->linear_solver); assert(direction.status == SubproblemStatus::OPTIMAL && "The primal-dual perturbed subproblem was not solved to optimality"); this->number_subproblems_solved++; @@ -209,13 +210,14 @@ namespace uno { } void PrimalDualInteriorPointSubproblem::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, - const Multipliers& current_multipliers) { + const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { // assemble, factorize and regularize the augmented matrix this->augmented_system.assemble_matrix(this->hessian_model->hessian, this->constraint_jacobian, problem.number_variables, problem.number_constraints); - this->augmented_system.factorize_matrix(problem.model, *this->linear_solver); + DEBUG << "Testing factorization with regularization factors (0, 0)\n"; + this->augmented_system.factorize_matrix(*this->linear_solver, warmstart_information); const double dual_regularization_parameter = std::pow(this->barrier_parameter(), this->parameters.regularization_exponent); - this->augmented_system.regularize_matrix(statistics, problem.model, *this->linear_solver, problem.number_variables, problem.number_constraints, - dual_regularization_parameter); + this->augmented_system.regularize_matrix(statistics, *this->linear_solver, problem.number_variables, problem.number_constraints, + dual_regularization_parameter, warmstart_information); // check the inertia [[maybe_unused]] auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = this->linear_solver->get_inertia(); diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp index bef8aefc..ee5401af 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.hpp @@ -79,7 +79,8 @@ namespace uno { const Vector& primal_direction, double tau); [[nodiscard]] static double dual_fraction_to_boundary(const OptimizationProblem& problem, const Multipliers& current_multipliers, Multipliers& direction_multipliers, double tau); - void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers); + void assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, const Multipliers& current_multipliers, + const WarmstartInformation& warmstart_information); void assemble_augmented_rhs(const OptimizationProblem& problem, const Multipliers& current_multipliers); void assemble_primal_dual_direction(const OptimizationProblem& problem, const Vector& current_primals, const Multipliers& current_multipliers, Vector& direction_primals, Multipliers& direction_multipliers); diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 69176875..37aa9bf8 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -10,8 +10,9 @@ #include "RectangularMatrix.hpp" #include "ingredients/hessian_models/UnstableRegularization.hpp" #include "model/Model.hpp" -#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" +#include "optimization/WarmstartInformation.hpp" #include "options/Options.hpp" +#include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "tools/Statistics.hpp" namespace uno { @@ -26,14 +27,13 @@ namespace uno { const Options& options); void assemble_matrix(const SymmetricMatrix& hessian, const RectangularMatrix& constraint_jacobian, size_t number_variables, size_t number_constraints); - void factorize_matrix(const Model& model, DirectSymmetricIndefiniteLinearSolver& linear_solver); - void regularize_matrix(Statistics& statistics, const Model& model, DirectSymmetricIndefiniteLinearSolver& linear_solver, - size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter); + void factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information); + void regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, + size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information); void solve(DirectSymmetricIndefiniteLinearSolver& linear_solver); // [[nodiscard]] T get_primal_regularization() const; protected: - size_t number_factorizations{0}; ElementType primal_regularization{0.}; ElementType dual_regularization{0.}; ElementType previous_primal_regularization{0.}; @@ -89,34 +89,36 @@ namespace uno { } template - void SymmetricIndefiniteLinearSystem::factorize_matrix(const Model& /*model*/, - DirectSymmetricIndefiniteLinearSolver& linear_solver) { - if (true || this->number_factorizations == 0) { + void SymmetricIndefiniteLinearSystem::factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver, + const WarmstartInformation& warmstart_information) { + if (warmstart_information.problem_structure_changed) { + DEBUG << "Performing symbolic analysis of the indefinite system\n"; linear_solver.do_symbolic_analysis(this->matrix); } + DEBUG << "Performing numerical factorization of the indefinite system\n"; linear_solver.do_numerical_factorization(this->matrix); - this->number_factorizations++; } + // the matrix has been factorized prior to calling this function template - void SymmetricIndefiniteLinearSystem::regularize_matrix(Statistics& statistics, const Model& model, + void SymmetricIndefiniteLinearSystem::regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, size_t size_primal_block, size_t size_dual_block, - ElementType dual_regularization_parameter) { + ElementType dual_regularization_parameter, const WarmstartInformation& warmstart_information) { DEBUG2 << "Original matrix\n" << this->matrix << '\n'; this->primal_regularization = ElementType(0.); this->dual_regularization = ElementType(0.); size_t number_attempts = 1; - DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n"; + DEBUG << "Number of attempts: " << number_attempts << "\n\n"; + + auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia(); + DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0)\n"; + DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n"; - if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) { - DEBUG << "Inertia is good\n"; + if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) { + DEBUG << "The inertia is correct\n"; statistics.set("regulariz", this->primal_regularization); return; } - auto [number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues] = linear_solver.get_inertia(); - DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), "; - DEBUG << "got (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n"; - DEBUG << "Number of attempts: " << number_attempts << "\n\n"; // set the constraint regularization coefficient if (linear_solver.matrix_is_singular()) { @@ -141,19 +143,20 @@ namespace uno { while (not good_inertia) { DEBUG << "Testing factorization with regularization factors (" << this->primal_regularization << ", " << this->dual_regularization << ")\n"; DEBUG2 << this->matrix << '\n'; - this->factorize_matrix(model, linear_solver); + this->factorize_matrix(linear_solver, warmstart_information); number_attempts++; + DEBUG << "Number of attempts: " << number_attempts << "\n"; + + std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia(); + DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0)\n"; + DEBUG << "Estimated inertia (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n"; - if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) { + if (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) { good_inertia = true; - DEBUG << "Factorization was a success\n"; + DEBUG << "The inertia is correct\n"; this->previous_primal_regularization = this->primal_regularization; } else { - std::tie(number_pos_eigenvalues, number_neg_eigenvalues, number_zero_eigenvalues) = linear_solver.get_inertia(); - DEBUG << "Expected inertia (" << size_primal_block << ", " << size_dual_block << ", 0), "; - DEBUG << "got (" << number_pos_eigenvalues << ", " << number_neg_eigenvalues << ", " << number_zero_eigenvalues << ")\n"; - DEBUG << "Number of attempts: " << number_attempts << "\n"; if (this->previous_primal_regularization == 0. || this->threshold_unsuccessful_attempts < number_attempts) { this->primal_regularization *= this->primal_regularization_fast_increase_factor; } @@ -188,4 +191,4 @@ namespace uno { */ } // namespace -#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H \ No newline at end of file +#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H diff --git a/uno/optimization/WarmstartInformation.cpp b/uno/optimization/WarmstartInformation.cpp index badabff9..1cefbb5a 100644 --- a/uno/optimization/WarmstartInformation.cpp +++ b/uno/optimization/WarmstartInformation.cpp @@ -10,7 +10,7 @@ namespace uno { std::cout << "Constraints: " << std::boolalpha << this->constraints_changed << '\n'; std::cout << "Constraint bounds: " << std::boolalpha << this->constraint_bounds_changed << '\n'; std::cout << "Variable bounds: " << std::boolalpha << this->variable_bounds_changed << '\n'; - std::cout << "Problem: " << std::boolalpha << this->problem_changed << '\n'; + std::cout << "Problem structure: " << std::boolalpha << this->problem_structure_changed << '\n'; } void WarmstartInformation::no_changes() { @@ -18,7 +18,7 @@ namespace uno { this->constraints_changed = false; this->constraint_bounds_changed = false; this->variable_bounds_changed = false; - this->problem_changed = false; + this->problem_structure_changed = false; } void WarmstartInformation::iterate_changed() { @@ -33,7 +33,7 @@ namespace uno { this->constraints_changed = true; this->constraint_bounds_changed = true; this->variable_bounds_changed = true; - this->problem_changed = true; + this->problem_structure_changed = true; } void WarmstartInformation::only_objective_changed() { @@ -41,6 +41,6 @@ namespace uno { this->constraints_changed = false; this->constraint_bounds_changed = false; this->variable_bounds_changed = false; - this->problem_changed = false; + this->problem_structure_changed = false; } } // namespace \ No newline at end of file diff --git a/uno/optimization/WarmstartInformation.hpp b/uno/optimization/WarmstartInformation.hpp index 0a6a9f01..0cb5fba1 100644 --- a/uno/optimization/WarmstartInformation.hpp +++ b/uno/optimization/WarmstartInformation.hpp @@ -10,7 +10,7 @@ namespace uno { bool constraints_changed{true}; bool constraint_bounds_changed{true}; bool variable_bounds_changed{true}; - bool problem_changed{true}; + bool problem_structure_changed{true}; void display() const; void no_changes(); diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index 58aeea7a..92eaf087 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/uno/solvers/BQPD/BQPDSolver.cpp @@ -165,8 +165,8 @@ namespace uno { BQPDMode BQPDSolver::determine_mode(const WarmstartInformation& warmstart_information) const { BQPDMode mode = (this->number_calls == 0) ? BQPDMode::ACTIVE_SET_EQUALITIES : BQPDMode::USER_DEFINED; - // if problem changed, use cold start - if (warmstart_information.problem_changed) { + // if problem structure changed, use cold start + if (warmstart_information.problem_structure_changed) { mode = BQPDMode::ACTIVE_SET_EQUALITIES; } // if only the variable bounds changed, reuse the active set estimate and the Jacobian information diff --git a/uno/solvers/MA57/MA57Solver.cpp b/uno/solvers/MA57/MA57Solver.cpp index 0688f66a..3b7b536f 100644 --- a/uno/solvers/MA57/MA57Solver.cpp +++ b/uno/solvers/MA57/MA57Solver.cpp @@ -165,7 +165,7 @@ namespace uno { // build the internal matrix representation this->row_indices.clear(); this->column_indices.clear(); - for (const auto [row_index, column_index, element]: matrix) { + for (const auto [row_index, column_index, _]: matrix) { this->row_indices.emplace_back(static_cast(row_index + this->fortran_shift)); this->column_indices.emplace_back(static_cast(column_index + this->fortran_shift)); } diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 6a26f257..791a0d78 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -29,25 +29,40 @@ namespace uno { this->mumps_structure.icntl[1] = -1; this->mumps_structure.icntl[2] = -1; this->mumps_structure.icntl[3] = 0; + //this->mumps_structure.icntl[5] = 0; // no scaling + //this->mumps_structure.icntl[7] = 0; // no scaling + this->mumps_structure.icntl[12] = 1; this->mumps_structure.icntl[23] = 1; // ICNTL(24) controls the detection of “null pivot rows” + + /* + // debug for MUMPS team + this->mumps_structure.icntl[1] = 6; // ICNTL(2)=6 + this->mumps_structure.icntl[2] = 6; // ICNTL(3)=6 + this->mumps_structure.icntl[3] = 6; // ICNTL(4)=2 + */ + + this->row_indices.reserve(number_nonzeros); + this->column_indices.reserve(number_nonzeros); } MUMPSSolver::~MUMPSSolver() { this->mumps_structure.job = MUMPSSolver::JOB_END; dmumps_c(&this->mumps_structure); } - + void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix& matrix) { - this->save_sparsity_to_local_format(matrix); + this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS; this->mumps_structure.n = static_cast(matrix.dimension()); this->mumps_structure.nnz = static_cast(matrix.number_nonzeros()); - this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS; this->mumps_structure.a = nullptr; - // connect the local COO matrix with the pointers in the structure + this->save_sparsity_to_local_format(matrix); + // connect the local sparsity with the pointers in the structure this->mumps_structure.irn = this->row_indices.data(); this->mumps_structure.jcn = this->column_indices.data(); + this->mumps_structure.a = nullptr; dmumps_c(&this->mumps_structure); + this->mumps_structure.icntl[7] = 8; // ICNTL(8) = 8: recompute scaling before factorization } void MUMPSSolver::do_numerical_factorization(const SymmetricMatrix& matrix) { @@ -93,8 +108,8 @@ namespace uno { this->row_indices.clear(); this->column_indices.clear(); for (const auto [row_index, column_index, _]: matrix) { - this->row_indices.emplace_back(row_index + this->fortran_shift); - this->column_indices.emplace_back(column_index + this->fortran_shift); + this->row_indices.emplace_back(static_cast(row_index + this->fortran_shift)); + this->column_indices.emplace_back(static_cast(column_index + this->fortran_shift)); } } } // namespace diff --git a/uno/solvers/MUMPS/MUMPSSolver.hpp b/uno/solvers/MUMPS/MUMPSSolver.hpp index da7ebd98..88aa143d 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.hpp +++ b/uno/solvers/MUMPS/MUMPSSolver.hpp @@ -27,6 +27,7 @@ namespace uno { protected: DMUMPS_STRUC_C mumps_structure{}; + // matrix sparsity std::vector row_indices{}; std::vector column_indices{};