From 1d786b627a79617b22217938ad9f198d128abc70 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Sun, 8 Dec 2024 00:42:31 +0100 Subject: [PATCH] Fixed the allocated memory for augmented system and linear solver in IPM class + updated the WarmstartInformation object in the regularization loop --- .../l1Relaxation.cpp | 4 ++-- .../l1Relaxation.hpp | 4 ++-- uno/ingredients/subproblems/Subproblem.hpp | 2 +- .../LPSubproblem.cpp | 2 +- .../LPSubproblem.hpp | 2 +- .../QPSubproblem.cpp | 2 +- .../QPSubproblem.hpp | 2 +- .../PrimalDualInteriorPointMethod.cpp | 15 ++++++++++----- .../PrimalDualInteriorPointMethod.hpp | 10 +++++----- .../SymmetricIndefiniteLinearSystem.hpp | 7 +++++-- uno/solvers/MA27/MA27Solver.cpp | 2 +- uno/solvers/MA57/MA57Solver.cpp | 6 +++--- 12 files changed, 33 insertions(+), 25 deletions(-) diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp index b49ed36b..4722d4b9 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.cpp @@ -136,7 +136,7 @@ namespace uno { } void l1Relaxation::solve_subproblem(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, - const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { + const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) { DEBUG << "Solving the subproblem with penalty parameter " << problem.get_objective_multiplier() << "\n\n"; // solve the subproblem @@ -148,7 +148,7 @@ namespace uno { } void l1Relaxation::solve_l1_relaxed_problem(Statistics& statistics, Iterate& current_iterate, Direction& direction, - double current_penalty_parameter, const WarmstartInformation& warmstart_information) { + double current_penalty_parameter, WarmstartInformation& warmstart_information) { this->l1_relaxed_problem.set_objective_multiplier(current_penalty_parameter); this->solve_subproblem(statistics, this->l1_relaxed_problem, current_iterate, current_iterate.multipliers, direction, warmstart_information); } diff --git a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp index e3635d0c..c10ae6a5 100644 --- a/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/l1Relaxation.hpp @@ -58,9 +58,9 @@ namespace uno { void solve_sequence_of_relaxed_subproblems(Statistics& statistics, Iterate& current_iterate, Direction& direction, WarmstartInformation& warmstart_information); void solve_l1_relaxed_problem(Statistics& statistics, Iterate& current_iterate, Direction& direction, double current_penalty_parameter, - const WarmstartInformation& warmstart_information); + WarmstartInformation& warmstart_information); void solve_subproblem(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information); + Direction& direction, WarmstartInformation& warmstart_information); // functions that decrease the penalty parameter to enforce particular conditions void decrease_parameter_aggressively(Iterate& current_iterate, const Direction& direction); diff --git a/uno/ingredients/subproblems/Subproblem.hpp b/uno/ingredients/subproblems/Subproblem.hpp index 8ff9ce92..6c355eff 100644 --- a/uno/ingredients/subproblems/Subproblem.hpp +++ b/uno/ingredients/subproblems/Subproblem.hpp @@ -37,7 +37,7 @@ namespace uno { virtual void initialize_statistics(Statistics& statistics, const Options& options) = 0; virtual void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) = 0; virtual void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) = 0; + Direction& direction, WarmstartInformation& warmstart_information) = 0; void set_trust_region_radius(double new_trust_region_radius); virtual void initialize_feasibility_problem(const l1RelaxedProblem& problem, Iterate& current_iterate) = 0; diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp index 7d3c9298..80b5d506 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.cpp @@ -36,7 +36,7 @@ namespace uno { } void LPSubproblem::solve(Statistics& /*statistics*/, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) { + Direction& direction, WarmstartInformation& warmstart_information) { // evaluate the functions at the current iterate this->evaluate_functions(problem, current_iterate, warmstart_information); diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp index b0670c0e..669aa0d0 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/LPSubproblem.hpp @@ -20,7 +20,7 @@ namespace uno { void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) override; + Direction& direction, WarmstartInformation& warmstart_information) override; private: // pointer to allow polymorphism diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp index 4c7e2511..6de22090 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp @@ -60,7 +60,7 @@ namespace uno { } void QPSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) { + Direction& direction, WarmstartInformation& warmstart_information) { // evaluate the functions at the current iterate this->evaluate_functions(statistics, problem, current_iterate, current_multipliers, warmstart_information); diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp index b0eea453..27f9f58a 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp +++ b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.hpp @@ -20,7 +20,7 @@ namespace uno { void initialize_statistics(Statistics& statistics, const Options& options) override; void generate_initial_iterate(const OptimizationProblem& problem, Iterate& initial_iterate) override; void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) override; + Direction& direction, WarmstartInformation& warmstart_information) override; protected: const bool use_regularization; diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp index e916aa6a..8e2c6149 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.cpp @@ -20,21 +20,26 @@ namespace uno { PrimalDualInteriorPointMethod::PrimalDualInteriorPointMethod(size_t number_variables, size_t number_constraints, size_t number_jacobian_nonzeros, size_t number_hessian_nonzeros, const Options& options): - Subproblem("exact", number_variables, number_hessian_nonzeros, false, options), + Subproblem("exact", number_variables, number_hessian_nonzeros + + number_variables /* diagonal barrier terms for bound constraints */ + + number_jacobian_nonzeros, /* Jacobian */ + false, options), objective_gradient(2 * number_variables), // original variables + barrier terms constraints(number_constraints), constraint_jacobian(number_constraints, number_variables), augmented_system(options.get_string("sparse_format"), number_variables + number_constraints, number_hessian_nonzeros + number_variables /* diagonal barrier terms for bound constraints */ - + number_jacobian_nonzeros /* Jacobian */, + + number_jacobian_nonzeros /* Jacobian */ + + number_variables + number_constraints, /* regularization */ true, /* use regularization */ options), linear_solver(SymmetricIndefiniteLinearSolverFactory::create(number_variables + number_constraints, number_hessian_nonzeros + number_variables + number_constraints /* regularization */ + 2 * number_variables /* diagonal barrier terms */ - + number_jacobian_nonzeros, /* Jacobian */ + + number_jacobian_nonzeros /* Jacobian */ + + number_variables + number_constraints, /* regularization */ options)), barrier_parameter_update_strategy(options), previous_barrier_parameter(options.get_double("barrier_initial_parameter")), @@ -161,7 +166,7 @@ namespace uno { } void PrimalDualInteriorPointMethod::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, - const Multipliers& current_multipliers, Direction& direction, const WarmstartInformation& warmstart_information) { + const Multipliers& current_multipliers, Direction& direction, WarmstartInformation& warmstart_information) { if (problem.has_inequality_constraints()) { throw std::runtime_error("The problem has inequality constraints. Create an instance of HomogeneousEqualityConstrainedModel"); } @@ -196,7 +201,7 @@ namespace uno { } void PrimalDualInteriorPointMethod::assemble_augmented_system(Statistics& statistics, const OptimizationProblem& problem, - const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) { + const Multipliers& current_multipliers, 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"; diff --git a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp index cbbf021b..7bde6fae 100644 --- a/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp +++ b/uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointMethod.hpp @@ -1,8 +1,8 @@ // Copyright (c) 2018-2024 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. -#ifndef UNO_INFEASIBLEINTERIORPOINTMETHOD_H -#define UNO_INFEASIBLEINTERIORPOINTMETHOD_H +#ifndef UNO_PRIMALDUALINTERIORPOINTMETHOD_H +#define UNO_PRIMALDUALINTERIORPOINTMETHOD_H #include "ingredients/subproblems/Subproblem.hpp" #include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp" @@ -39,7 +39,7 @@ namespace uno { void exit_feasibility_problem(const OptimizationProblem& problem, Iterate& trial_iterate) override; void solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate, const Multipliers& current_multipliers, - Direction& direction, const WarmstartInformation& warmstart_information) override; + Direction& direction, WarmstartInformation& warmstart_information) override; void set_auxiliary_measure(const Model& model, Iterate& iterate) override; [[nodiscard]] double compute_predicted_auxiliary_reduction_model(const Model& model, const Iterate& current_iterate, @@ -81,7 +81,7 @@ namespace uno { [[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, - const WarmstartInformation& warmstart_information); + 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); @@ -91,4 +91,4 @@ namespace uno { }; } // namespace -#endif // UNO_INFEASIBLEINTERIORPOINTMETHOD_H \ No newline at end of file +#endif // UNO_PRIMALDUALINTERIORPOINTMETHOD_H \ No newline at end of file diff --git a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp index 25f9accc..68781ce5 100644 --- a/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp +++ b/uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp @@ -30,7 +30,7 @@ namespace uno { size_t number_variables, size_t number_constraints); 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); + size_t size_primal_block, size_t size_dual_block, ElementType dual_regularization_parameter, WarmstartInformation& warmstart_information); void solve(DirectSymmetricIndefiniteLinearSolver& linear_solver); // [[nodiscard]] T get_primal_regularization() const; @@ -104,7 +104,10 @@ namespace uno { template void SymmetricIndefiniteLinearSystem::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) { + ElementType dual_regularization_parameter, WarmstartInformation& warmstart_information) { + // the linear system was already analyzed by the linear solver + warmstart_information.problem_structure_changed = false; + DEBUG2 << "Original matrix\n" << this->matrix << '\n'; this->primal_regularization = ElementType(0.); this->dual_regularization = ElementType(0.); diff --git a/uno/solvers/MA27/MA27Solver.cpp b/uno/solvers/MA27/MA27Solver.cpp index 41041f0f..57277b02 100644 --- a/uno/solvers/MA27/MA27Solver.cpp +++ b/uno/solvers/MA27/MA27Solver.cpp @@ -108,7 +108,6 @@ namespace uno { RANK_DEFICIENT, // Matrix is rank deficient. In this case, a decomposition will still have been produced which will enable the subsequent solution of consistent equations (MA27B/BD entry only). INFO(2) will be set to the rank of the matrix. Note that this warning will overwrite an INFO(1)=1 or INFO(1)=2 warning. }; - MA27Solver::MA27Solver(size_t max_dimension, size_t max_number_nonzeros): DirectSymmetricIndefiniteLinearSolver(max_dimension), n(static_cast(max_dimension)), nnz(static_cast(max_number_nonzeros)), @@ -125,6 +124,7 @@ namespace uno { icntl[eICNTL::LDIAG] = 0; } + void MA27Solver::do_symbolic_analysis(const SymmetricMatrix& matrix) { assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size"); assert(matrix.number_nonzeros() <= irn.capacity() && diff --git a/uno/solvers/MA57/MA57Solver.cpp b/uno/solvers/MA57/MA57Solver.cpp index 3b7b536f..805a85f3 100644 --- a/uno/solvers/MA57/MA57Solver.cpp +++ b/uno/solvers/MA57/MA57Solver.cpp @@ -73,9 +73,9 @@ namespace uno { /* out */ this->info.data(), /* out */ this->rinfo.data()); - assert(0 <= info[0] && "MA57: the symbolic analysis failed"); - if (0 < info[0]) { - WARNING << "MA57 has issued a warning: info(1) = " << info[0] << '\n'; + assert(0 <= this->info[0] && "MA57: the symbolic analysis failed"); + if (0 < this->info[0]) { + WARNING << "MA57 has issued a warning: info(1) = " << this->info[0] << '\n'; } // get LFACT and LIFACT and resize FACT and IFACT (no effect if resized to <= size)