From 2462129f7e81eaad4967ebb31a82f717c16ef91d Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 19 Nov 2024 01:12:55 +0100 Subject: [PATCH 01/11] Removed fixed_hessian_sparsity flag in Model --- .../PrimalDualInteriorPointSubproblem.cpp | 15 ++++++------ .../PrimalDualInteriorPointSubproblem.hpp | 3 ++- .../SymmetricIndefiniteLinearSystem.hpp | 23 ++++++++++--------- 3 files changed, 22 insertions(+), 19 deletions(-) diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index 06a07436..5896e9d8 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" @@ -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,13 @@ 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); + 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..2bbfbd57 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,9 +27,9 @@ 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; @@ -89,9 +90,9 @@ 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_changed) { linear_solver.do_symbolic_analysis(this->matrix); } linear_solver.do_numerical_factorization(this->matrix); @@ -99,9 +100,9 @@ namespace uno { } 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.); @@ -141,7 +142,7 @@ 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++; if (not linear_solver.matrix_is_singular() && linear_solver.number_negative_eigenvalues() == size_dual_block) { @@ -188,4 +189,4 @@ namespace uno { */ } // namespace -#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H \ No newline at end of file +#endif // UNO_SYMMETRICINDEFINITELINEARSYSTEM_H From 11af9bbfc1f1a6f0d1f79c94bfc9c6e3f132c977 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 19 Nov 2024 12:18:20 +0100 Subject: [PATCH 02/11] Minor refactoring --- .../PrimalDualInteriorPointSubproblem.cpp | 1 + .../SymmetricIndefiniteLinearSystem.hpp | 28 ++++++++++--------- 2 files changed, 16 insertions(+), 13 deletions(-) diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index 5896e9d8..8a2cb158 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp @@ -213,6 +213,7 @@ namespace uno { 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); + 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, *this->linear_solver, problem.number_variables, problem.number_constraints, diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 2bbfbd57..4a3d9faf 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -99,6 +99,7 @@ namespace uno { this->number_factorizations++; } + // the matrix has been factorized prior to calling this function template void SymmetricIndefiniteLinearSystem::regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver& linear_solver, size_t size_primal_block, size_t size_dual_block, @@ -107,17 +108,17 @@ namespace uno { 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), "; + 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()) { @@ -144,17 +145,18 @@ namespace uno { DEBUG2 << this->matrix << '\n'; 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), "; + 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; } From efc3fd560835c36fdbf82187675fcdc039546b3b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 19 Nov 2024 13:01:35 +0100 Subject: [PATCH 03/11] Pass warmstart object to factorize_matrix() (without using it for the moment) + minor refactoring --- uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp | 6 +++--- uno/optimization/WarmstartInformation.cpp | 8 ++++---- uno/optimization/WarmstartInformation.hpp | 2 +- uno/solvers/BQPD/BQPDSolver.cpp | 4 ++-- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 4a3d9faf..07d530b9 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -34,7 +34,6 @@ namespace uno { // [[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.}; @@ -92,11 +91,12 @@ namespace uno { template void SymmetricIndefiniteLinearSystem::factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information) { - if (warmstart_information.problem_changed) { + 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 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 From d74fa6f7f6f6485fd09adb55175e0e17acb65016 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 19 Nov 2024 14:16:56 +0100 Subject: [PATCH 04/11] Trying out other MUMPS parameters --- .../SymmetricIndefiniteLinearSystem.hpp | 4 ++-- uno/solvers/MUMPS/MUMPSSolver.cpp | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 07d530b9..37aa9bf8 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -111,7 +111,7 @@ namespace uno { 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), "; + 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 (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) { @@ -148,7 +148,7 @@ namespace uno { 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), "; + 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 (number_pos_eigenvalues == size_primal_block && number_neg_eigenvalues == size_dual_block && number_zero_eigenvalues == 0) { diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 6a26f257..ee83c5a2 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -26,9 +26,12 @@ namespace uno { dmumps_c(&this->mumps_structure); // control parameters this->mumps_structure.icntl[0] = -1; - this->mumps_structure.icntl[1] = -1; - this->mumps_structure.icntl[2] = -1; + this->mumps_structure.icntl[1] = 0; + this->mumps_structure.icntl[2] = 0; this->mumps_structure.icntl[3] = 0; + this->mumps_structure.icntl[6] = 7; // pivot order + this->mumps_structure.icntl[7] = 6; // scaling + this->mumps_structure.icntl[9] = 0; // no iterative refinement this->mumps_structure.icntl[12] = 1; this->mumps_structure.icntl[23] = 1; // ICNTL(24) controls the detection of “null pivot rows” } @@ -37,12 +40,12 @@ namespace uno { this->mumps_structure.job = MUMPSSolver::JOB_END; dmumps_c(&this->mumps_structure); } - + void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix& matrix) { + this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS; this->save_sparsity_to_local_format(matrix); 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->mumps_structure.irn = this->row_indices.data(); From fadf598ab57fa08e2d2919ba651374f10b7f029f Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 19 Nov 2024 22:13:00 +0100 Subject: [PATCH 05/11] Restored previous MUMPS settings --- uno/solvers/MUMPS/MUMPSSolver.cpp | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index ee83c5a2..05b58277 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -26,12 +26,9 @@ namespace uno { dmumps_c(&this->mumps_structure); // control parameters this->mumps_structure.icntl[0] = -1; - this->mumps_structure.icntl[1] = 0; - this->mumps_structure.icntl[2] = 0; + this->mumps_structure.icntl[1] = -1; + this->mumps_structure.icntl[2] = -1; this->mumps_structure.icntl[3] = 0; - this->mumps_structure.icntl[6] = 7; // pivot order - this->mumps_structure.icntl[7] = 6; // scaling - this->mumps_structure.icntl[9] = 0; // no iterative refinement this->mumps_structure.icntl[12] = 1; this->mumps_structure.icntl[23] = 1; // ICNTL(24) controls the detection of “null pivot rows” } From de1c6d68dc10e9030db619a83c24d304d47cab93 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 25 Nov 2024 10:42:43 +0100 Subject: [PATCH 06/11] MUMPSSolver: replaced member COO_matrix with two vectors --- uno/solvers/MA57/MA57Solver.cpp | 2 +- uno/solvers/MUMPS/MUMPSSolver.cpp | 12 ++++++++---- uno/solvers/MUMPS/MUMPSSolver.hpp | 1 + 3 files changed, 10 insertions(+), 5 deletions(-) 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 05b58277..64efe0d3 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -29,8 +29,12 @@ namespace uno { this->mumps_structure.icntl[1] = -1; this->mumps_structure.icntl[2] = -1; this->mumps_structure.icntl[3] = 0; + this->mumps_structure.icntl[12] = 1; this->mumps_structure.icntl[23] = 1; // ICNTL(24) controls the detection of “null pivot rows” + + this->row_indices.reserve(number_nonzeros); + this->column_indices.reserve(number_nonzeros); } MUMPSSolver::~MUMPSSolver() { @@ -40,11 +44,11 @@ namespace uno { void MUMPSSolver::do_symbolic_analysis(const SymmetricMatrix& matrix) { this->mumps_structure.job = MUMPSSolver::JOB_ANALYSIS; - this->save_sparsity_to_local_format(matrix); this->mumps_structure.n = static_cast(matrix.dimension()); this->mumps_structure.nnz = static_cast(matrix.number_nonzeros()); 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(); dmumps_c(&this->mumps_structure); @@ -93,8 +97,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{}; From 87ecfdd517642fd0d4ec153891adab17c066a58b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 25 Nov 2024 10:45:39 +0100 Subject: [PATCH 07/11] MUMPSSolver: disable scaling during analysis --- uno/solvers/MUMPS/MUMPSSolver.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 64efe0d3..b17a267e 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -29,6 +29,8 @@ 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” @@ -51,6 +53,7 @@ namespace uno { // 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); } From 97acb41dbf4d5ecf306574163955f96460cb9864 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 25 Nov 2024 10:47:40 +0100 Subject: [PATCH 08/11] Relaxed threshold for unbounded problems --- .github/julia/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 From 5f08f5937d36be3d5563d529246fa31f20183fa4 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Mon, 25 Nov 2024 10:48:03 +0100 Subject: [PATCH 09/11] SymmetricIndefiniteLinearSystem: perform the analysis only when the problem has changed --- .../PrimalDualInteriorPointSubproblem.cpp | 40 +++++++++---------- uno/solvers/MUMPS/MUMPSSolver.cpp | 4 +- 2 files changed, 22 insertions(+), 22 deletions(-) diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp index 8a2cb158..1a58b3c7 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp @@ -124,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 @@ -175,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, diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index b17a267e..713f473a 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -29,8 +29,8 @@ 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[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” From 30c7304367c84360af9d1e67c057cc525dbc4025 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Tue, 3 Dec 2024 23:06:22 +0100 Subject: [PATCH 10/11] Temporarily restore analysis at each iteration --- uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp | 2 +- uno/solvers/MUMPS/MUMPSSolver.cpp | 7 +++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 37aa9bf8..b53d072c 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -91,7 +91,7 @@ namespace uno { template void SymmetricIndefiniteLinearSystem::factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information) { - if (warmstart_information.problem_structure_changed) { + if (true) { DEBUG << "Performing symbolic analysis of the indefinite system\n"; linear_solver.do_symbolic_analysis(this->matrix); } diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 713f473a..62d6482f 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -35,6 +35,13 @@ namespace uno { 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); } From 5670b8eb6c8c2b0f51283f8b7ca82f097bc0292b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 4 Dec 2024 23:46:35 +0100 Subject: [PATCH 11/11] Perform analysis only when the problem structure changes + set a setting in MUMPS to recompute scaling before each factorization --- uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp | 2 +- uno/solvers/MUMPS/MUMPSSolver.cpp | 1 + 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index b53d072c..37aa9bf8 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -91,7 +91,7 @@ namespace uno { template void SymmetricIndefiniteLinearSystem::factorize_matrix(DirectSymmetricIndefiniteLinearSolver& linear_solver, const WarmstartInformation& warmstart_information) { - if (true) { + if (warmstart_information.problem_structure_changed) { DEBUG << "Performing symbolic analysis of the indefinite system\n"; linear_solver.do_symbolic_analysis(this->matrix); } diff --git a/uno/solvers/MUMPS/MUMPSSolver.cpp b/uno/solvers/MUMPS/MUMPSSolver.cpp index 62d6482f..791a0d78 100644 --- a/uno/solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/solvers/MUMPS/MUMPSSolver.cpp @@ -62,6 +62,7 @@ namespace uno { 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) {