Skip to content

Commit

Permalink
Minor refactoring
Browse files Browse the repository at this point in the history
  • Loading branch information
@cvanaret
cvanaret committed Nov 19, 2024
1 parent 306489a commit 8f098dd
Show file tree
Hide file tree
Showing 2 changed files with 23 additions and 21 deletions.
5 changes: 3 additions & 2 deletions uno/ingredients/subproblems/interior_point_methods/PrimalDualInteriorPointSubproblem.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,9 +212,10 @@ namespace uno {
const Multipliers& current_multipliers) {
// 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);
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,
this->augmented_system.regularize_matrix(statistics, *this->linear_solver, problem.number_variables, problem.number_constraints,
dual_regularization_parameter);

// check the inertia
Expand Down
39 changes: 20 additions & 19 deletions uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,8 +26,8 @@ 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(const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
void regularize_matrix(Statistics& statistics, const Model& model, DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver,
void factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
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);
void solve(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver);
// [[nodiscard]] T get_primal_regularization() const;
Expand Down Expand Up @@ -89,8 +89,7 @@ namespace uno {
}

template <typename ElementType>
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(const Model& model,
DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver) {
void SymmetricIndefiniteLinearSystem<ElementType>::factorize_matrix(DirectSymmetricIndefiniteLinearSolver<size_t, ElementType>& linear_solver) {
// compute the symbolic factorization only when:
// the problem has a non-constant augmented system (ie is not an LP or a QP) or it is the first factorization
if (true || this->number_factorizations == 0) {
Expand All @@ -100,25 +99,26 @@ namespace uno {
this->number_factorizations++;
}

// the matrix has been factorized prior to calling this function
template <typename ElementType>
void SymmetricIndefiniteLinearSystem<ElementType>::regularize_matrix(Statistics& statistics, const Model& model,
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) {
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), ";
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()) {
Expand All @@ -143,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);
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;
}
Expand Down

0 comments on commit 8f098dd

Please sign in to comment.