Fixed the allocated memory for augmented system and linear solver in …

…IPM class + updated the WarmstartInformation object in the regularization loop

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/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";

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<double>& current_primals, const Multipliers& current_multipliers,
             Vector<double>& direction_primals, Multipliers& direction_multipliers);
@@ -91,4 +91,4 @@ namespace uno {
    };
 } // namespace
 
-#endif // UNO_INFEASIBLEINTERIORPOINTMETHOD_H
+#endif // UNO_PRIMALDUALINTERIORPOINTMETHOD_H

uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp

@@ -30,7 +30,7 @@ namespace uno {
             size_t number_variables, size_t number_constraints);
       void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver, const 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;
 
@@ -104,7 +104,10 @@ 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) {
+      // 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.);

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<size_t, double>(max_dimension),
          n(static_cast<int>(max_dimension)), nnz(static_cast<int>(max_number_nonzeros)),
@@ -125,6 +124,7 @@ namespace uno {
       icntl[eICNTL::LDIAG] = 0;
    }
 
+
    void MA27Solver::do_symbolic_analysis(const SymmetricMatrix<size_t, double>& matrix) {
       assert(matrix.dimension() <= iw1.capacity() && "MA27Solver: the dimension of the matrix is larger than the preallocated size");
       assert(matrix.number_nonzeros() <= irn.capacity() &&

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)