Merge pull request #129 from cvanaret/indefinite_system_allocation

Fixed memory allocation in IPM + removed useless symbolic analyses

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);
    }

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);

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;

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);
 

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

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);
 

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;

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp

@@ -3,7 +3,6 @@
 
 #include <cmath>
 #include "PrimalDualInteriorPointSubproblem.hpp"
-#include "ingredients/hessian_models/HessianModelFactory.hpp"
 #include "linear_algebra/SparseStorageFactory.hpp"
 #include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp"
 #include "solvers/DirectSymmetricIndefiniteLinearSolver.hpp"
@@ -26,14 +25,16 @@ namespace uno {
          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")),
@@ -178,7 +179,7 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointSubproblem::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");
       }
@@ -210,7 +211,7 @@ namespace uno {
    }
 
    void PrimalDualInteriorPointSubproblem::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";

uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.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_INFEASIBLEINTERIORPOINTSUBPROBLEM_H
-#define UNO_INFEASIBLEINTERIORPOINTSUBPROBLEM_H
+#ifndef UNO_PRIMALDUALINTERIORPOINTSUBPROBLEM_H
+#define UNO_PRIMALDUALINTERIORPOINTSUBPROBLEM_H
 
 #include "ingredients/subproblems/Subproblem.hpp"
 #include "linear_algebra/SymmetricIndefiniteLinearSystem.hpp"
@@ -38,7 +38,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,
@@ -80,7 +80,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<double>& current_primals, const Multipliers& current_multipliers,
             Vector<double>& direction_primals, Multipliers& direction_multipliers);
@@ -90,4 +90,4 @@ namespace uno {
    };
 } // namespace
 
-#endif // UNO_INFEASIBLEINTERIORPOINTSUBPROBLEM_H
+#endif // UNO_PRIMALDUALINTERIORPOINTSUBPROBLEM_H

uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp

@@ -27,9 +27,9 @@ namespace uno {
             const Options& options);
       void assemble_matrix(const SymmetricMatrix<size_t, double>& hessian, const RectangularMatrix<double>& constraint_jacobian,
             size_t number_variables, size_t number_constraints);
-      void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, const WarmstartInformation& warmstart_information);
+      void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, WarmstartInformation& warmstart_information);
       void regularize_matrix(Statistics& statistics, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& 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<size_t, ElementType>& linear_solver);
       // [[nodiscard]] T get_primal_regularization() const;
 
@@ -90,10 +90,11 @@ namespace uno {
 
    template <typename ElementType>
    void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
-         const WarmstartInformation& warmstart_information) {
+         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);
+         warmstart_information.problem_structure_changed = false;
       }
       DEBUG << "Performing numerical factorization of the indefinite system\n";
       linear_solver.do_numerical_factorization(this->matrix);
@@ -103,7 +104,7 @@ namespace uno {
    template <typename ElementType>
    void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics,
          DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& 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) {
       DEBUG2 << "Original matrix\n" << this->matrix << '\n';
       this->primal_regularization = ElementType(0.);
       this->dual_regularization = ElementType(0.);