From 657e35a26068d73b7749ee9f2400bd93351b32a6 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 17:07:01 +0200 Subject: [PATCH 01/13] Replaced raw pointer with View (renamed VectorView) in compute_hessian_vector_product --- interfaces/AMPL/AMPLModel.cpp | 6 +- interfaces/AMPL/AMPLModel.hpp | 4 +- interfaces/C/Uno_C_API.cpp | 7 ++- interfaces/Python/cpp_classes/PythonModel.cpp | 10 ++-- interfaces/Python/cpp_classes/PythonModel.hpp | 4 +- .../ConstraintRelaxationStrategy.cpp | 2 +- .../FeasibilityRestoration.cpp | 2 +- .../NoRelaxation.cpp | 2 +- .../relaxed_problems/l1RelaxedProblem.cpp | 6 +- .../relaxed_problems/l1RelaxedProblem.hpp | 4 +- .../hessian_models/ExactHessian.cpp | 4 +- .../hessian_models/ExactHessian.hpp | 4 +- .../hessian_models/HessianModel.hpp | 5 +- .../hessian_models/IdentityHessian.cpp | 4 +- .../hessian_models/IdentityHessian.hpp | 4 +- .../hessian_models/ZeroHessian.cpp | 4 +- .../hessian_models/ZeroHessian.hpp | 4 +- .../direct/DirectQuasiNewtonHessian.hpp | 2 +- .../quasi_newton/direct/LBFGSHessian.cpp | 12 ++-- .../quasi_newton/direct/LBFGSHessian.hpp | 6 +- .../quasi_newton/direct/LSR1Hessian.cpp | 8 +-- .../quasi_newton/direct/LSR1Hessian.hpp | 6 +- .../inverse/InverseLBFGSHessian.cpp | 6 +- .../inverse/InverseLBFGSHessian.hpp | 4 +- .../PrimalDualInteriorPointProblem.cpp | 6 +- .../PrimalDualInteriorPointProblem.hpp | 4 +- uno/ingredients/subproblem/Subproblem.cpp | 4 +- uno/ingredients/subproblem/Subproblem.hpp | 4 +- .../subproblem_solvers/BQPD/BQPDSolver.cpp | 14 +++-- .../subproblem_solvers/BQPD/BQPDWorkspace.cpp | 6 +- .../subproblem_solvers/SubproblemSolver.cpp | 2 +- uno/linear_algebra/BLASMatrix.hpp | 2 +- uno/linear_algebra/DenseMatrix.hpp | 10 ++-- uno/linear_algebra/Vector.hpp | 6 +- .../{VectorView.hpp => View.hpp} | 59 ++++++++++--------- uno/model/BoundRelaxedModel.hpp | 4 +- uno/model/FixedBoundsConstraintsModel.cpp | 6 +- uno/model/FixedBoundsConstraintsModel.hpp | 4 +- .../HomogeneousEqualityConstrainedModel.cpp | 6 +- .../HomogeneousEqualityConstrainedModel.hpp | 4 +- uno/model/Model.cpp | 2 +- uno/model/Model.hpp | 5 +- uno/optimization/Direction.cpp | 2 +- uno/optimization/OptimizationProblem.cpp | 6 +- uno/optimization/OptimizationProblem.hpp | 7 ++- uno/optimization/Result.cpp | 2 +- 46 files changed, 147 insertions(+), 138 deletions(-) rename uno/linear_algebra/{VectorView.hpp => View.hpp} (85%) diff --git a/interfaces/AMPL/AMPLModel.cpp b/interfaces/AMPL/AMPLModel.cpp index 3dc9ed9cf..2767311f1 100644 --- a/interfaces/AMPL/AMPLModel.cpp +++ b/interfaces/AMPL/AMPLModel.cpp @@ -225,13 +225,13 @@ namespace uno { throw std::runtime_error("AMPLModel::compute_jacobian_transposed_vector_product not implemented"); } - void AMPLModel::compute_hessian_vector_product(const double* /*x*/, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const { + void AMPLModel::compute_hessian_vector_product(View /*x*/, View vector, double objective_multiplier, + const Vector& multipliers, View result) const { // scale by the objective sign objective_multiplier *= this->optimization_sense; // compute the Hessian-vector product - (this->asl->p.Hvcomp)(this->asl, result, const_cast(vector), -1, &objective_multiplier, + (this->asl->p.Hvcomp)(this->asl, result.data(), const_cast(vector.data()), -1, &objective_multiplier, const_cast(multipliers.data())); } diff --git a/interfaces/AMPL/AMPLModel.hpp b/interfaces/AMPL/AMPLModel.hpp index 52750a190..331b3b591 100644 --- a/interfaces/AMPL/AMPLModel.hpp +++ b/interfaces/AMPL/AMPLModel.hpp @@ -57,8 +57,8 @@ namespace uno { // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override; [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; diff --git a/interfaces/C/Uno_C_API.cpp b/interfaces/C/Uno_C_API.cpp index d2d956f60..c08df82a3 100644 --- a/interfaces/C/Uno_C_API.cpp +++ b/interfaces/C/Uno_C_API.cpp @@ -210,8 +210,8 @@ class UnoModel: public Model { } } - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override { + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override { if (this->user_model.lagrangian_hessian_operator != nullptr) { objective_multiplier *= this->optimization_sense; // if the model has a different sign convention for the Lagrangian than Uno, flip the signs of the multipliers @@ -219,7 +219,8 @@ class UnoModel: public Model { const_cast&>(multipliers).scale(-1.); } const uno_int return_code = this->user_model.lagrangian_hessian_operator(this->user_model.number_variables, - this->user_model.number_constraints, x, true, objective_multiplier, multipliers.data(), vector, result, this->user_model.user_data); + this->user_model.number_constraints, x.data(), true, objective_multiplier, multipliers.data(), vector.data(), + result.data(), this->user_model.user_data); // flip the signs of the multipliers back if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { const_cast&>(multipliers).scale(-1.); diff --git a/interfaces/Python/cpp_classes/PythonModel.cpp b/interfaces/Python/cpp_classes/PythonModel.cpp index 8b6ac969d..5231b4dd8 100644 --- a/interfaces/Python/cpp_classes/PythonModel.cpp +++ b/interfaces/Python/cpp_classes/PythonModel.cpp @@ -96,7 +96,7 @@ namespace uno { } } - void PythonModel::compute_jacobian_sparsity(uno_int * row_indices, uno_int * column_indices, uno_int row_offset, + void PythonModel::compute_jacobian_sparsity(uno_int* row_indices, uno_int* column_indices, uno_int row_offset, uno_int column_offset, uno_int solver_indexing, MatrixOrder /*matrix_format*/) const { // copy the indices of the user sparsity patterns to the Uno vectors for (size_t index: Range(static_cast(this->user_model.number_jacobian_nonzeros))) { @@ -219,8 +219,8 @@ namespace uno { } } - void PythonModel::compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const { + void PythonModel::compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const { if (this->user_model.lagrangian_hessian_operator.has_value()) { objective_multiplier *= this->optimization_sense; // if the model has a different sign convention for the Lagrangian than Uno, flip the signs of the multipliers @@ -229,8 +229,8 @@ namespace uno { } const auto x_py = to_const_array(x, this->number_variables); const auto multipliers_py = to_const_array(multipliers.data(), this->number_constraints); - const auto vector_py = to_const_array(vector, this->number_variables); - auto result_py = to_array(result, this->number_variables); + const auto vector_py = to_const_array(vector.data(), this->number_variables); + auto result_py = to_array(result.data(), this->number_variables); // evaluate Hessian-vector product try { diff --git a/interfaces/Python/cpp_classes/PythonModel.hpp b/interfaces/Python/cpp_classes/PythonModel.hpp index 04b0c150c..bf1083609 100644 --- a/interfaces/Python/cpp_classes/PythonModel.hpp +++ b/interfaces/Python/cpp_classes/PythonModel.hpp @@ -47,8 +47,8 @@ namespace uno { // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override; // purely functions [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; diff --git a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp index 799261cf3..27880360c 100644 --- a/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/ConstraintRelaxationStrategy.cpp @@ -4,7 +4,7 @@ #include "ConstraintRelaxationStrategy.hpp" #include "ingredients/globalization_strategies/GlobalizationStrategy.hpp" #include "ingredients/subproblem/Subproblem.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/EvaluationCache.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 7b8f1df36..7b9dfc8fb 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -14,7 +14,7 @@ #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategyFactory.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "ingredients/subproblem_solvers/SubproblemSolverFactory.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/EvaluationCache.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp index 39d35fe17..8fcea7c00 100644 --- a/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/NoRelaxation.cpp @@ -11,7 +11,7 @@ #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategyFactory.hpp" #include "ingredients/subproblem/Subproblem.hpp" #include "ingredients/subproblem_solvers/SubproblemSolverFactory.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Direction.hpp" #include "optimization/EvaluationCache.hpp" #include "optimization/Iterate.hpp" diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp index fbd3a7f7a..133eab763 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp @@ -3,7 +3,7 @@ #include "l1RelaxedProblem.hpp" #include "ingredients/hessian_models/HessianModel.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "model/Model.hpp" #include "optimization/Evaluations.hpp" #include "optimization/Iterate.hpp" @@ -263,8 +263,8 @@ namespace uno { } } - void l1RelaxedProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const { + void l1RelaxedProblem::compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, + const Multipliers& multipliers, View result) const { hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result); // proximal contribution diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp index fd0894d86..b39a6e742 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp @@ -46,8 +46,8 @@ namespace uno { // linear operators void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; void compute_jacobian_transposed_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; - void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const override; + void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, + const Multipliers& multipliers, View result) const override; [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; diff --git a/uno/ingredients/hessian_models/ExactHessian.cpp b/uno/ingredients/hessian_models/ExactHessian.cpp index e7cd84a68..414450fbf 100644 --- a/uno/ingredients/hessian_models/ExactHessian.cpp +++ b/uno/ingredients/hessian_models/ExactHessian.cpp @@ -46,8 +46,8 @@ namespace uno { ++this->evaluation_count; } - void ExactHessian::compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& constraint_multipliers, double* result) { + void ExactHessian::compute_hessian_vector_product(View x, View vector, + double objective_multiplier, const Vector& constraint_multipliers, View result) { this->model.compute_hessian_vector_product(x, vector, objective_multiplier, constraint_multipliers, result); ++this->evaluation_count; } diff --git a/uno/ingredients/hessian_models/ExactHessian.hpp b/uno/ingredients/hessian_models/ExactHessian.hpp index 0f72558af..2cf845134 100644 --- a/uno/ingredients/hessian_models/ExactHessian.hpp +++ b/uno/ingredients/hessian_models/ExactHessian.hpp @@ -27,8 +27,8 @@ namespace uno { EvaluationCache& evaluation_cache) override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& constraint_multipliers, View result) override; protected: const Model& model; diff --git a/uno/ingredients/hessian_models/HessianModel.hpp b/uno/ingredients/hessian_models/HessianModel.hpp index 258faa0db..ef3f31d24 100644 --- a/uno/ingredients/hessian_models/HessianModel.hpp +++ b/uno/ingredients/hessian_models/HessianModel.hpp @@ -8,6 +8,7 @@ #include #include #include "../interfaces/C/uno_int.h" +#include "linear_algebra/View.hpp" namespace uno { // forward declarations @@ -37,8 +38,8 @@ namespace uno { EvaluationCache& evaluation_cache) = 0; virtual void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) = 0; - virtual void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) = 0; + virtual void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& constraint_multipliers, View result) = 0; }; } // namespace diff --git a/uno/ingredients/hessian_models/IdentityHessian.cpp b/uno/ingredients/hessian_models/IdentityHessian.cpp index bb7b9544c..565b22212 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.cpp +++ b/uno/ingredients/hessian_models/IdentityHessian.cpp @@ -52,8 +52,8 @@ namespace uno { } } - void IdentityHessian::compute_hessian_vector_product(const double* /*x*/, const double* vector, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result) { + void IdentityHessian::compute_hessian_vector_product(View /*x*/, View vector, + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, View result) { for (size_t variable_index: Range(this->number_variables)) { result[variable_index] = vector[variable_index]; } diff --git a/uno/ingredients/hessian_models/IdentityHessian.hpp b/uno/ingredients/hessian_models/IdentityHessian.hpp index 84ad28742..1417bab27 100644 --- a/uno/ingredients/hessian_models/IdentityHessian.hpp +++ b/uno/ingredients/hessian_models/IdentityHessian.hpp @@ -23,8 +23,8 @@ namespace uno { EvaluationCache& evaluation_cache) override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& constraint_multipliers, View result) override; protected: const size_t number_variables; diff --git a/uno/ingredients/hessian_models/ZeroHessian.cpp b/uno/ingredients/hessian_models/ZeroHessian.cpp index 0f6442826..6667146a6 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.cpp +++ b/uno/ingredients/hessian_models/ZeroHessian.cpp @@ -44,8 +44,8 @@ namespace uno { // do nothing } - void ZeroHessian::compute_hessian_vector_product(const double* /*x*/, const double* /*vector*/, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* result) { + void ZeroHessian::compute_hessian_vector_product(View /*x*/, View /*vector*/, + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, View result) { for (size_t variable_index: Range(this->number_variables)) { result[variable_index] = 0.; } diff --git a/uno/ingredients/hessian_models/ZeroHessian.hpp b/uno/ingredients/hessian_models/ZeroHessian.hpp index e081b2bf8..28fc51232 100644 --- a/uno/ingredients/hessian_models/ZeroHessian.hpp +++ b/uno/ingredients/hessian_models/ZeroHessian.hpp @@ -24,8 +24,8 @@ namespace uno { EvaluationCache& evaluation_cache) override; void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& constraint_multipliers, View result) override; protected: const size_t number_variables; diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/DirectQuasiNewtonHessian.hpp b/uno/ingredients/hessian_models/quasi_newton/direct/DirectQuasiNewtonHessian.hpp index ffbbd1499..923ddb5ae 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/DirectQuasiNewtonHessian.hpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/DirectQuasiNewtonHessian.hpp @@ -24,7 +24,7 @@ namespace uno { // functions that can be called by WoodburyEQPSolver [[nodiscard]] virtual size_t get_correction_rank() const = 0; - [[nodiscard]] virtual VectorView get_correction_column(size_t column_index) const = 0; + [[nodiscard]] virtual View get_correction_column(size_t column_index) const = 0; [[nodiscard]] virtual double get_correction_column_scaling(size_t column_index) const = 0; }; } // namespace diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp index b8c0caeaa..78c759d9b 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp @@ -6,7 +6,7 @@ #include "model/Model.hpp" #include "options/Options.hpp" #include "linear_algebra/BLAS.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" #include "symbolic/Inverse.hpp" #include "symbolic/Multiplication.hpp" @@ -85,8 +85,8 @@ namespace uno { // Hessian-vector product where the Hessian approximation is Bk = B0 - U Uᵀ + V Vᵀ and B0 = δ I // Bk v = (B0 - U Uᵀ + V Vᵀ) v = δ v - U (Uᵀ v) + V (Vᵀ v) - void LBFGSHessian::compute_hessian_vector_product(const double* /*x*/, const double* vector, - double objective_multiplier, const Vector& /*constraint_multipliers*/, double* result) { + void LBFGSHessian::compute_hessian_vector_product(View /*x*/, View vector, + double objective_multiplier, const Vector& /*constraint_multipliers*/, View result) { if (objective_multiplier != this->fixed_objective_multiplier) { throw std::runtime_error("The quasi-Newton Hessian model was initialized with a different objective multiplier"); } @@ -109,14 +109,14 @@ namespace uno { const auto current_U_column = this->U.column(column_index); const double U_coefficient = -dot(current_U_column, vector); // minus sign for U // result += coefficient * current_column - blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result); + blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result.data()); } // work on each column of V (Vᵀ v) for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_V_column = this->V.column(column_index); const double V_coefficient = dot(current_V_column, vector); // plus sign for V // result += coefficient * current_column - blas1::add(this->model.number_variables, V_coefficient, current_V_column.data(), result); + blas1::add(this->model.number_variables, V_coefficient, current_V_column.data(), result.data()); } } @@ -125,7 +125,7 @@ namespace uno { } // get a column of (U V) - VectorView LBFGSHessian::get_correction_column(size_t column_index) const { + View LBFGSHessian::get_correction_column(size_t column_index) const { if (column_index < this->number_entries_in_memory) { return this->U.column(column_index); } diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.hpp b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.hpp index 15ac1dc0d..327189b9f 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.hpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.hpp @@ -29,12 +29,12 @@ namespace uno { void notify_trial_iterate(Statistics& statistics, const Iterate& current_iterate, const Iterate& trial_iterate, EvaluationCache& evaluation_cache) override; - void compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, + double objective_multiplier, const Vector& constraint_multipliers, View result) override; // functions that can be called by WoodburyEQPSolver [[nodiscard]] size_t get_correction_rank() const override; - [[nodiscard]] VectorView get_correction_column(size_t column_index) const override; + [[nodiscard]] View get_correction_column(size_t column_index) const override; [[nodiscard]] double get_correction_column_scaling(size_t column_index) const override; protected: diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp index a4f2952dd..e5d8d5b2b 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp @@ -53,8 +53,8 @@ namespace uno { // Hessian-vector product where the Hessian approximation is Bk = B0 + U P⁻¹ Uᵀ and B0 = δ I // Bk v = (B0 + U P⁻¹ Uᵀ) v = δ v + U (P⁻¹ (Uᵀ v)) - void LSR1Hessian::compute_hessian_vector_product(const double* /*x*/, const double* vector, - double objective_multiplier, const Vector& /*constraint_multipliers*/, double* result) { + void LSR1Hessian::compute_hessian_vector_product(View /*x*/, View vector, + double objective_multiplier, const Vector& /*constraint_multipliers*/, View result) { if (objective_multiplier != this->fixed_objective_multiplier) { throw std::runtime_error("The L-SR1 Hessian model was initialized with a different objective multiplier"); } @@ -78,7 +78,7 @@ namespace uno { double U_coefficient = dot(current_U_column, vector) / this->get_correction_column_scaling(column_index); assert(!std::isnan(U_coefficient)); // result += coefficient * column(P⁻¹) * current_column - blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result); + blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result.data()); } } @@ -86,7 +86,7 @@ namespace uno { return this->number_entries_in_memory; } - VectorView LSR1Hessian::get_correction_column(size_t column_index) const { + View LSR1Hessian::get_correction_column(size_t column_index) const { return this->U.column(column_index); } diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.hpp b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.hpp index 047638e8d..fffa1dd04 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.hpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.hpp @@ -24,12 +24,12 @@ namespace uno { void notify_trial_iterate(Statistics& statistics, const Iterate& current_iterate, const Iterate& trial_iterate, EvaluationCache& evaluation_cache) override; - void compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, + double objective_multiplier, const Vector& constraint_multipliers, View result) override; // functions that can be called by WoodburyEQPSolver [[nodiscard]] size_t get_correction_rank() const override; - [[nodiscard]] VectorView get_correction_column(size_t column_index) const override; + [[nodiscard]] View get_correction_column(size_t column_index) const override; [[nodiscard]] double get_correction_column_scaling(size_t column_index) const override; protected: diff --git a/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.cpp b/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.cpp index 5f4e5f2e5..fd05bff5a 100644 --- a/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.cpp @@ -6,7 +6,7 @@ #include "model/Model.hpp" #include "options/Options.hpp" #include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" #include "symbolic/Inverse.hpp" #include "symbolic/Multiplication.hpp" @@ -80,8 +80,8 @@ namespace uno { throw std::runtime_error("InverseLBFGSHessian::evaluate_hessian should not be called."); } - void InverseLBFGSHessian::compute_hessian_vector_product(const double* /*x*/, const double* /*vector*/, - double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, double* /*result*/) { + void InverseLBFGSHessian::compute_hessian_vector_product(View /*x*/, View /*vector*/, + double /*objective_multiplier*/, const Vector& /*constraint_multipliers*/, View /*result*/) { throw std::runtime_error("InverseLBFGSHessian::compute_hessian_vector_product should not be called."); } diff --git a/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.hpp b/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.hpp index 520591e3f..1305a14a8 100644 --- a/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.hpp +++ b/uno/ingredients/hessian_models/quasi_newton/inverse/InverseLBFGSHessian.hpp @@ -28,8 +28,8 @@ namespace uno { void evaluate_hessian(Statistics& statistics, const Vector& primal_variables, double objective_multiplier, const Vector& constraint_multipliers, double* hessian_values) override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& constraint_multipliers, double* result) override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& constraint_multipliers, View result) override; // function that can be called by NewtonSolver void compute_inverse_hessian_vector_product(const double* x, const double* vector, double* result); diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp index 7bee3ba5f..2b4409506 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp @@ -5,7 +5,7 @@ #include "../InteriorPointParameters.hpp" #include "ingredients/hessian_models/HessianModel.hpp" #include "linear_algebra/SparseVector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/Evaluations.hpp" @@ -242,8 +242,8 @@ namespace uno { this->inner.compute_jacobian_transposed_vector_product(vector, result, evaluations); } - void PrimalDualInteriorPointProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, - const double* vector, const Multipliers& multipliers, double* result) const { + void PrimalDualInteriorPointProblem::compute_hessian_vector_product(HessianModel& hessian_model, View x, + View vector, const Multipliers& multipliers, View result) const { // original Lagrangian Hessian this->inner.compute_hessian_vector_product(hessian_model, x, vector, multipliers, result); diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp index 683c6aee9..42e0dc905 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp @@ -46,8 +46,8 @@ namespace uno { // linear operators void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; void compute_jacobian_transposed_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; - void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const override; + void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, + const Multipliers& multipliers, View result) const override; [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index e1292d55b..efd128768 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -6,8 +6,6 @@ #include "ingredients/inertia_correction_strategies/InertiaCorrectionStrategy.hpp" #include "ingredients/subproblem_solvers/DirectSymmetricIndefiniteLinearSolver.hpp" #include "ingredients/subproblem_solvers/SolverWorkspace.hpp" -#include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/Evaluations.hpp" @@ -99,7 +97,7 @@ namespace uno { } } - void Subproblem::compute_hessian_vector_product(const double* x, const double* vector, double* result) const { + void Subproblem::compute_hessian_vector_product(View x, View vector, View result) const { // unregularized Hessian-vector product this->problem.compute_hessian_vector_product(this->hessian_model, x, vector, this->current_iterate.multipliers, result); diff --git a/uno/ingredients/subproblem/Subproblem.hpp b/uno/ingredients/subproblem/Subproblem.hpp index 94b103800..7a1a01244 100644 --- a/uno/ingredients/subproblem/Subproblem.hpp +++ b/uno/ingredients/subproblem/Subproblem.hpp @@ -8,7 +8,7 @@ #include "ingredients/globalization_strategies/ProgressMeasures.hpp" #include "linear_algebra/MatrixOrder.hpp" #include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/OptimizationProblem.hpp" #include "symbolic/IntegerRange.hpp" @@ -41,7 +41,7 @@ namespace uno { // regularized Hessian void evaluate_lagrangian_hessian(Statistics& statistics, double* hessian_values) const; void regularize_lagrangian_hessian(Statistics& statistics, double* hessian_values) const; - void compute_hessian_vector_product(const double* x, const double* vector, double* result) const; + void compute_hessian_vector_product(View x, View vector, View result) const; // augmented system void regularize_augmented_matrix(Statistics& statistics, double* augmented_matrix_values, diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index 4cfb046c3..2fbb13248 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -7,7 +7,7 @@ #include "ingredients/subproblem/Subproblem.hpp" #include "linear_algebra/Indexing.hpp" #include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Direction.hpp" #include "optimization/Iterate.hpp" #include "optimization/WarmstartInformation.hpp" @@ -20,7 +20,7 @@ #define hessian_vector_product FC_GLOBAL(gdotx, GDOTX) extern "C" { - void hessian_vector_product([[maybe_unused]] int *dimension, const double vector[], const double ws[], + void hessian_vector_product([[maybe_unused]] int *int_dimension, const double vector[], const double ws[], const int lws[], double result[]); // fortran common block used in bqpd/bqpd.f @@ -329,10 +329,11 @@ namespace uno { } } // namespace -void hessian_vector_product(int* dimension, const double vector[], const double /*ws*/[], const int lws[], double result[]) { - assert(dimension != nullptr && "BQPDSolver::hessian_vector_product: the dimension n passed by pointer is NULL"); +void hessian_vector_product(int* int_dimension, const double vector[], const double /*ws*/[], const int lws[], double result[]) { + assert(int_dimension != nullptr && *int_dimension >= 0); + const size_t dimension = static_cast(*int_dimension); - for (size_t i = 0; i < static_cast(*dimension); i++) { + for (size_t i = 0; i < dimension; i++) { result[i] = 0.; } @@ -372,6 +373,7 @@ void hessian_vector_product(int* dimension, const double vector[], const double } // otherwise, try to perform a Hessian-vector product if possible else if (subproblem->has_hessian_operator()) { - subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.data(), vector, result); + subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.view(), + uno::view(vector, dimension), uno::view(result, dimension)); } } \ No newline at end of file diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDWorkspace.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDWorkspace.cpp index ef0cf0b5f..68d0ab60f 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDWorkspace.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDWorkspace.cpp @@ -7,7 +7,7 @@ #include "ingredients/subproblem/Subproblem.hpp" #include "linear_algebra/Indexing.hpp" #include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" #include "optimization/WarmstartInformation.hpp" @@ -43,8 +43,8 @@ namespace uno { if (subproblem.has_hessian_operator()) { // linear operator // TODO compute the quadratic form directly without temporary result // compute Hv - subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.data(), vector.data(), - this->hessian_vector_product.data()); + subproblem.compute_hessian_vector_product(subproblem.current_iterate.primals.view(), vector.view(), + this->hessian_vector_product.view()); // compute the dot product return dot(view(vector, 0, subproblem.number_variables), this->hessian_vector_product); } diff --git a/uno/ingredients/subproblem_solvers/SubproblemSolver.cpp b/uno/ingredients/subproblem_solvers/SubproblemSolver.cpp index 6791153bc..376d90e4e 100644 --- a/uno/ingredients/subproblem_solvers/SubproblemSolver.cpp +++ b/uno/ingredients/subproblem_solvers/SubproblemSolver.cpp @@ -3,7 +3,7 @@ #include "SubproblemSolver.hpp" #include "ingredients/subproblem/Subproblem.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" #include "optimization/Multipliers.hpp" diff --git a/uno/linear_algebra/BLASMatrix.hpp b/uno/linear_algebra/BLASMatrix.hpp index e02d1f335..099912bec 100644 --- a/uno/linear_algebra/BLASMatrix.hpp +++ b/uno/linear_algebra/BLASMatrix.hpp @@ -7,7 +7,7 @@ #include "Vector.hpp" #include "linear_algebra/BLAS.hpp" #include "linear_algebra/LAPACK.hpp" -#include "VectorView.hpp" +#include "View.hpp" #include "symbolic/Inverse.hpp" #include "symbolic/Multiplication.hpp" #include "symbolic/Sum.hpp" diff --git a/uno/linear_algebra/DenseMatrix.hpp b/uno/linear_algebra/DenseMatrix.hpp index 093a34403..73a055892 100644 --- a/uno/linear_algebra/DenseMatrix.hpp +++ b/uno/linear_algebra/DenseMatrix.hpp @@ -7,7 +7,7 @@ #include #include #include "linear_algebra/BLASMatrix.hpp" -#include "VectorView.hpp" +#include "View.hpp" #include "symbolic/Range.hpp" namespace uno { @@ -62,8 +62,8 @@ namespace uno { [[nodiscard]] T& entry(size_t row_index, size_t column_index); [[nodiscard]] const T& entry(size_t row_index, size_t column_index) const; // vector view - [[nodiscard]] VectorView column(size_t column_index); - [[nodiscard]] VectorView column(size_t column_index) const; + [[nodiscard]] View column(size_t column_index); + [[nodiscard]] View column(size_t column_index) const; [[nodiscard]] Submatrix submatrix(size_t number_rows, size_t number_columns); [[nodiscard]] T* data() override; @@ -92,12 +92,12 @@ namespace uno { } template - VectorView DenseMatrix::column(size_t column_index) { + View DenseMatrix::column(size_t column_index) { return {this->matrix.data() + column_index * this->number_rows, this->number_rows}; } template - VectorView DenseMatrix::column(size_t column_index) const { + View DenseMatrix::column(size_t column_index) const { return {this->matrix.data() + column_index * this->number_rows, this->number_rows}; } diff --git a/uno/linear_algebra/Vector.hpp b/uno/linear_algebra/Vector.hpp index 53f3b2df7..62a0f549e 100644 --- a/uno/linear_algebra/Vector.hpp +++ b/uno/linear_algebra/Vector.hpp @@ -7,7 +7,7 @@ #include #include #include -#include "VectorView.hpp" +#include "View.hpp" #include "symbolic/Range.hpp" namespace uno { @@ -149,11 +149,11 @@ namespace uno { this->view().scale(factor); } - VectorView view() const noexcept { + View view() const noexcept { return uno::view(this->data(), this->size()); } - VectorView view() noexcept { + View view() noexcept { return uno::view(this->data(), this->size()); } diff --git a/uno/linear_algebra/VectorView.hpp b/uno/linear_algebra/View.hpp similarity index 85% rename from uno/linear_algebra/VectorView.hpp rename to uno/linear_algebra/View.hpp index a32c9d8e4..cf6f01eec 100644 --- a/uno/linear_algebra/VectorView.hpp +++ b/uno/linear_algebra/View.hpp @@ -1,8 +1,8 @@ // Copyright (c) 2026 Charlie Vanaret // Licensed under the MIT license. See LICENSE file in the project directory for details. -#ifndef UNO_VECTORVIEW_H -#define UNO_VECTORVIEW_H +#ifndef UNO_VIEW_H +#define UNO_VIEW_H #include #include "BLAS.hpp" @@ -17,7 +17,7 @@ namespace uno { // constant contiguous array in memory on which BLAS can be called template - class VectorView { + class View { protected: T* pointer; const size_t view_size; @@ -25,11 +25,16 @@ namespace uno { public: using value_type = std::remove_const_t; - VectorView(T* pointer, size_t size): pointer(pointer), view_size(size) { } - ~VectorView() = default; - VectorView(const VectorView& other) = default; - VectorView(VectorView&& other) = default; - VectorView& operator=(const VectorView& other) { + View(T* pointer, size_t size): pointer(pointer), view_size(size) { } + ~View() = default; + View(const View& other) = default; + + // conversion View -> View + template >> + View(const View& other) noexcept: pointer(other.data()), view_size(other.size()) { } + + View(View&& other) = default; + View& operator=(const View& other) { if (&other != this) { blas1::copy(this->size(), other.data(), this->data()); } @@ -62,7 +67,7 @@ namespace uno { // specialized operation y = x, when the other vector has the member function data() template - VectorView& operator=(const Vector& other) { + View& operator=(const Vector& other) { if (other.size() != this->size()) { throw std::invalid_argument("Dimension mismatch between x and y"); } @@ -72,7 +77,7 @@ namespace uno { // generic operation y = expression template - VectorView& operator=(const Expression& expression) { + View& operator=(const Expression& expression) { // static_assert(std::is_same_v); if (expression.size() != this->size()) { throw std::invalid_argument("Dimension mismatch between expression and y"); @@ -85,7 +90,7 @@ namespace uno { // specialized operation y = a * x (note: no BLAS operation available) template - VectorView& operator=(ScalarMultiple&& expression) { + View& operator=(ScalarMultiple&& expression) { const auto& a = expression.get_factor(); const auto& x = expression.get_expression(); if (x.size() != this->size()) { @@ -100,7 +105,7 @@ namespace uno { // specialized operation y = x - z // note: no BLAS operation available template - VectorView& operator=(Subtraction&& expression) { + View& operator=(Subtraction&& expression) { const auto& x = expression.get_left(); const auto& z = expression.get_right(); if (x.size() != z.size()) { @@ -116,7 +121,7 @@ namespace uno { // y = Ax (or y = A x + 0 * y) template - VectorView& operator=(Multiplication&& expression) { + View& operator=(Multiplication&& expression) { const auto& x = expression.get_right(); if (x.data() == this->data()) { throw std::invalid_argument("BLASVector::operator= cannot be called with x == y"); @@ -135,7 +140,7 @@ namespace uno { // specialized operation y += a * x template - VectorView& operator+=(ScalarMultiple&& expression) { + View& operator+=(ScalarMultiple&& expression) { const auto& a = expression.get_factor(); const auto& x = expression.get_expression(); if (x.size() != this->size()) { @@ -147,7 +152,7 @@ namespace uno { // generic operation y += x template - VectorView& operator+=(const Vector& other) { + View& operator+=(const Vector& other) { if (other.size() != this->size()) { throw std::invalid_argument("Dimension mismatch between x and y"); } @@ -157,7 +162,7 @@ namespace uno { // generic operation y -= x template - VectorView& operator-=(const Vector& other) { + View& operator-=(const Vector& other) { if (other.size() != this->size()) { throw std::invalid_argument("Dimension mismatch between x and y"); } @@ -167,7 +172,7 @@ namespace uno { // x := U⁻ᵀ y with U upper triangular (solve Uᵀ x := y) template - VectorView& operator=(Multiplication>>, Vector>&& expression) { + View& operator=(Multiplication>>, Vector>&& expression) { const auto& U = expression.get_left().get_matrix().get_matrix().get_matrix(); const auto& y = expression.get_right(); if (this->size() != U.number_columns) { @@ -182,7 +187,7 @@ namespace uno { // x := U⁻¹ y with U upper triangular (solve U x := y) template - VectorView& operator=(Multiplication>, Vector>&& expression) { + View& operator=(Multiplication>, Vector>&& expression) { const auto& U = expression.get_left().get_matrix().get_matrix(); const auto& y = expression.get_right(); if (this->size() != U.number_columns) { @@ -197,7 +202,7 @@ namespace uno { // y := A^T x template - VectorView& operator=(Multiplication, Vector>&& expression) { + View& operator=(Multiplication, Vector>&& expression) { const auto& A = expression.get_left().get_matrix(); const auto& x = expression.get_right(); if (A.number_rows != x.size()) { @@ -213,7 +218,7 @@ namespace uno { // y += Ax (or y = A x + y) template - VectorView& operator+=(Multiplication&& expression) { + View& operator+=(Multiplication&& expression) { const auto& A = expression.get_left(); const auto& x = expression.get_right(); if (A.number_columns != x.size()) { @@ -229,7 +234,7 @@ namespace uno { // y -= Ax (or y = -A x + y) template - VectorView& operator-=(Multiplication&& expression) { + View& operator-=(Multiplication&& expression) { const auto& A = expression.get_left(); const auto& x = expression.get_right(); if (A.number_columns != x.size()) { @@ -262,12 +267,12 @@ namespace uno { if (start > end || end > container.size()) { throw std::out_of_range("invalid vector view"); } - return VectorView{container.data() + start, end - start}; + return View{container.data() + start, end - start}; } template auto view(T* pointer, size_t size) { - return VectorView{pointer, size}; + return View{pointer, size}; } template @@ -275,11 +280,11 @@ namespace uno { if (start > end) { throw std::out_of_range("invalid vector view"); } - return VectorView{pointer + start, end - start}; + return View{pointer + start, end - start}; } template - std::ostream& operator<<(std::ostream& stream, const VectorView& view) { + std::ostream& operator<<(std::ostream& stream, const View& view) { for (size_t index: Range(view.size())) { stream << view[index] << " "; } @@ -294,9 +299,9 @@ namespace uno { return blas1::dot(x.size(), x.data(), y.data()); } - inline double dot(const VectorView& x, const double* y) { + inline double dot(const View& x, const double* y) { return blas1::dot(x.size(), x.data(), y); } } // namespace -#endif //UNO_VECTORVIEW_H \ No newline at end of file +#endif //UNO_VIEW_H \ No newline at end of file diff --git a/uno/model/BoundRelaxedModel.hpp b/uno/model/BoundRelaxedModel.hpp index a83560653..8ce46fee7 100644 --- a/uno/model/BoundRelaxedModel.hpp +++ b/uno/model/BoundRelaxedModel.hpp @@ -77,8 +77,8 @@ namespace uno { this->model.compute_jacobian_transposed_vector_product(x, vector, result); } - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override { + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override { this->model.compute_hessian_vector_product(x, vector, objective_multiplier, multipliers, result); } diff --git a/uno/model/FixedBoundsConstraintsModel.cpp b/uno/model/FixedBoundsConstraintsModel.cpp index d9f25943b..8e407776b 100644 --- a/uno/model/FixedBoundsConstraintsModel.cpp +++ b/uno/model/FixedBoundsConstraintsModel.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "FixedBoundsConstraintsModel.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" namespace uno { @@ -133,8 +133,8 @@ namespace uno { } } - void FixedBoundsConstraintsModel::compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& multipliers, double* result) const { + void FixedBoundsConstraintsModel::compute_hessian_vector_product(View x, View vector, + double objective_multiplier, const Vector& multipliers, View result) const { this->model.compute_hessian_vector_product(x, vector, objective_multiplier, multipliers, result); } diff --git a/uno/model/FixedBoundsConstraintsModel.hpp b/uno/model/FixedBoundsConstraintsModel.hpp index b956ad5bf..acb136c13 100644 --- a/uno/model/FixedBoundsConstraintsModel.hpp +++ b/uno/model/FixedBoundsConstraintsModel.hpp @@ -48,8 +48,8 @@ namespace uno { // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override; [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; diff --git a/uno/model/HomogeneousEqualityConstrainedModel.cpp b/uno/model/HomogeneousEqualityConstrainedModel.cpp index 8d9d7cb8c..ecefd7ccf 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.cpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "HomogeneousEqualityConstrainedModel.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Iterate.hpp" #include "symbolic/Range.hpp" @@ -134,8 +134,8 @@ namespace uno { } } - void HomogeneousEqualityConstrainedModel::compute_hessian_vector_product(const double* x, const double* vector, - double objective_multiplier, const Vector& multipliers, double* result) const { + void HomogeneousEqualityConstrainedModel::compute_hessian_vector_product(View x, View vector, + double objective_multiplier, const Vector& multipliers, View result) const { this->model.compute_hessian_vector_product(x, vector, objective_multiplier, multipliers, result); } diff --git a/uno/model/HomogeneousEqualityConstrainedModel.hpp b/uno/model/HomogeneousEqualityConstrainedModel.hpp index 7bd31f81d..2cad620f8 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.hpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.hpp @@ -45,8 +45,8 @@ namespace uno { // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; - void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const override; + void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const override; [[nodiscard]] const std::vector& get_variables_lower_bounds() const override; [[nodiscard]] const std::vector& get_variables_upper_bounds() const override; diff --git a/uno/model/Model.cpp b/uno/model/Model.cpp index 54ede19e1..e2efaa739 100644 --- a/uno/model/Model.cpp +++ b/uno/model/Model.cpp @@ -6,7 +6,7 @@ #include #include "Model.hpp" #include "linear_algebra/Vector.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/Evaluations.hpp" #include "optimization/Iterate.hpp" #include "symbolic/ScalarMultiple.hpp" diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp index 04a9494d4..cee1f2c58 100644 --- a/uno/model/Model.hpp +++ b/uno/model/Model.hpp @@ -8,6 +8,7 @@ #include #include "linear_algebra/MatrixOrder.hpp" #include "linear_algebra/Norm.hpp" +#include "linear_algebra/View.hpp" #include "optimization/ProblemType.hpp" #include "symbolic/VectorExpression.hpp" #include "../interfaces/C/uno_int.h" @@ -74,8 +75,8 @@ namespace uno { // here we use pointers, since the vector and the result may be provided by a low-level subproblem solver virtual void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const = 0; virtual void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const = 0; - virtual void compute_hessian_vector_product(const double* x, const double* vector, double objective_multiplier, - const Vector& multipliers, double* result) const = 0; + virtual void compute_hessian_vector_product(View x, View vector, double objective_multiplier, + const Vector& multipliers, View result) const = 0; // purely virtual functions [[nodiscard]] virtual const std::vector& get_variables_lower_bounds() const = 0; diff --git a/uno/optimization/Direction.cpp b/uno/optimization/Direction.cpp index 745d7a92c..ff5ad1914 100644 --- a/uno/optimization/Direction.cpp +++ b/uno/optimization/Direction.cpp @@ -2,7 +2,7 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include "Direction.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "tools/Logger.hpp" namespace uno { diff --git a/uno/optimization/OptimizationProblem.cpp b/uno/optimization/OptimizationProblem.cpp index 452078e3b..1b18c663b 100644 --- a/uno/optimization/OptimizationProblem.cpp +++ b/uno/optimization/OptimizationProblem.cpp @@ -5,7 +5,7 @@ #include "ingredients/hessian_models/HessianModel.hpp" #include "ingredients/inequality_handling_methods/InequalityHandlingMethod.hpp" #include "linear_algebra/MatrixOrder.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "model/Model.hpp" #include "optimization/Direction.hpp" #include "optimization/Evaluations.hpp" @@ -114,8 +114,8 @@ namespace uno { evaluations.compute_jacobian_transposed_vector_product(this->model, vector, result); } - void OptimizationProblem::compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const { + void OptimizationProblem::compute_hessian_vector_product(HessianModel& hessian_model, View x, + View vector, const Multipliers& multipliers, View result) const { hessian_model.compute_hessian_vector_product(x, vector, this->get_objective_multiplier(), multipliers.constraints, result); } diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp index 564cd40cf..0c2e1c85c 100644 --- a/uno/optimization/OptimizationProblem.hpp +++ b/uno/optimization/OptimizationProblem.hpp @@ -9,6 +9,7 @@ #include "ingredients/inertia_correction_strategies/Inertia.hpp" #include "linear_algebra/MatrixOrder.hpp" #include "linear_algebra/Norm.hpp" +#include "linear_algebra/View.hpp" #include "optimization/SolutionStatus.hpp" #include "../interfaces/C/uno_int.h" #include "model/Model.hpp" @@ -35,7 +36,7 @@ namespace uno { explicit OptimizationProblem(const Model& model); OptimizationProblem(const Model& model, size_t number_variables, size_t number_constraints); virtual ~OptimizationProblem() = default; - virtual std::unique_ptr clone() const; + [[nodiscard]] virtual std::unique_ptr clone() const; const Model& model; const size_t number_variables; /*!< Number of variables */ @@ -70,8 +71,8 @@ namespace uno { virtual void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const; virtual void compute_jacobian_transposed_vector_product(const double* vector, double* result, const Evaluations& evaluations) const; - virtual void compute_hessian_vector_product(HessianModel& hessian_model, const double* x, const double* vector, - const Multipliers& multipliers, double* result) const; + virtual void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, + const Multipliers& multipliers, View result) const; [[nodiscard]] size_t get_number_original_variables() const; [[nodiscard]] virtual const std::vector& get_variables_lower_bounds() const; diff --git a/uno/optimization/Result.cpp b/uno/optimization/Result.cpp index 7fb8cf130..c2fe516c4 100644 --- a/uno/optimization/Result.cpp +++ b/uno/optimization/Result.cpp @@ -4,7 +4,7 @@ #include #include "Result.hpp" #include "SolutionStatus.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "tools/Logger.hpp" namespace uno { From 78a93096a085cbb22425531dc4aebd0d97987bff Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 17:27:09 +0200 Subject: [PATCH 02/13] Fixes --- interfaces/Python/cpp_classes/PythonModel.cpp | 2 +- unotest/unit_tests/{VectorViewTests.cpp => ViewTests.cpp} | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) rename unotest/unit_tests/{VectorViewTests.cpp => ViewTests.cpp} (85%) diff --git a/interfaces/Python/cpp_classes/PythonModel.cpp b/interfaces/Python/cpp_classes/PythonModel.cpp index 5231b4dd8..7e9911aae 100644 --- a/interfaces/Python/cpp_classes/PythonModel.cpp +++ b/interfaces/Python/cpp_classes/PythonModel.cpp @@ -5,7 +5,7 @@ #include #include #include "PythonModel.hpp" -#include "linear_algebra/VectorView.hpp" +#include "linear_algebra/View.hpp" #include "optimization/EvaluationErrors.hpp" #include "symbolic/Concatenation.hpp" #include "Uno.hpp" diff --git a/unotest/unit_tests/VectorViewTests.cpp b/unotest/unit_tests/ViewTests.cpp similarity index 85% rename from unotest/unit_tests/VectorViewTests.cpp rename to unotest/unit_tests/ViewTests.cpp index d74504940..673a70ee5 100644 --- a/unotest/unit_tests/VectorViewTests.cpp +++ b/unotest/unit_tests/ViewTests.cpp @@ -2,11 +2,12 @@ // Licensed under the MIT license. See LICENSE file in the project directory for details. #include -#include "linear_algebra/VectorView.hpp" +#include +#include "linear_algebra/View.hpp" using namespace uno; -TEST(VectorView, Size) { +TEST(View, Size) { const std::vector x{1., 2., 3., 100., 200., 300.}; // vector view const size_t start_index = 2; From 78639ecd27bf95e136f413cdcec277cc535c0714 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 17:33:41 +0200 Subject: [PATCH 03/13] Fix --- interfaces/Python/cpp_classes/PythonModel.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/interfaces/Python/cpp_classes/PythonModel.cpp b/interfaces/Python/cpp_classes/PythonModel.cpp index 7e9911aae..c6d5ac6a4 100644 --- a/interfaces/Python/cpp_classes/PythonModel.cpp +++ b/interfaces/Python/cpp_classes/PythonModel.cpp @@ -227,7 +227,7 @@ namespace uno { if (this->user_model.lagrangian_sign_convention == UNO_MULTIPLIER_POSITIVE) { const_cast&>(multipliers).scale(-1.); } - const auto x_py = to_const_array(x, this->number_variables); + const auto x_py = to_const_array(x.data(), this->number_variables); const auto multipliers_py = to_const_array(multipliers.data(), this->number_constraints); const auto vector_py = to_const_array(vector.data(), this->number_variables); auto result_py = to_array(result.data(), this->number_variables); From 7697191ff6d2a7606f6196f00773c192281974d0 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 19:17:09 +0200 Subject: [PATCH 04/13] Set logger to INFO --- .github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl | 2 +- .github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl b/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl index a56f56ca6..83b92dabc 100644 --- a/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl +++ b/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl @@ -22,7 +22,7 @@ function Optimizer(options) return AmplNLWriter.Optimizer(Uno_jll.amplexe, options) end -Optimizer_Uno_filtersqp() = Optimizer(["logger=SILENT", "preset=filtersqp", "QP_solver=BQPD", "globalization_mechanism=LS", "hessian_model=LBFGS", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) +Optimizer_Uno_filtersqp() = Optimizer(["logger=INFO", "preset=filtersqp", "QP_solver=BQPD", "globalization_mechanism=LS", "hessian_model=LBFGS", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) # This testset runs https://github.com/jump-dev/MINLPTests.jl @testset "MINLPTests" begin diff --git a/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl b/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl index cb337fda4..f51bc2239 100644 --- a/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl +++ b/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl @@ -22,7 +22,7 @@ function Optimizer(options) return AmplNLWriter.Optimizer(Uno_jll.amplexe, options) end -Optimizer_Uno_filtersqp() = Optimizer(["logger=SILENT", "preset=filtersqp", "QP_solver=BQPD", "hessian_model=LSR1", "quasi_newton_memory_size=20", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) +Optimizer_Uno_filtersqp() = Optimizer(["logger=INFO", "preset=filtersqp", "QP_solver=BQPD", "hessian_model=LSR1", "quasi_newton_memory_size=20", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) # This testset runs https://github.com/jump-dev/MINLPTests.jl @testset "MINLPTests" begin From 04b68ebeb6907d538b717ded9f9ecbad75fe8ad5 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 19:29:21 +0200 Subject: [PATCH 05/13] Create view of the right size --- .github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl | 2 +- .github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl | 2 +- uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl b/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl index 83b92dabc..a56f56ca6 100644 --- a/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl +++ b/.github/julia/runtests_uno_filtersqp_bqpd_ls_lbfgs.jl @@ -22,7 +22,7 @@ function Optimizer(options) return AmplNLWriter.Optimizer(Uno_jll.amplexe, options) end -Optimizer_Uno_filtersqp() = Optimizer(["logger=INFO", "preset=filtersqp", "QP_solver=BQPD", "globalization_mechanism=LS", "hessian_model=LBFGS", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) +Optimizer_Uno_filtersqp() = Optimizer(["logger=SILENT", "preset=filtersqp", "QP_solver=BQPD", "globalization_mechanism=LS", "hessian_model=LBFGS", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) # This testset runs https://github.com/jump-dev/MINLPTests.jl @testset "MINLPTests" begin diff --git a/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl b/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl index f51bc2239..cb337fda4 100644 --- a/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl +++ b/.github/julia/runtests_uno_filtersqp_bqpd_lsr1.jl @@ -22,7 +22,7 @@ function Optimizer(options) return AmplNLWriter.Optimizer(Uno_jll.amplexe, options) end -Optimizer_Uno_filtersqp() = Optimizer(["logger=INFO", "preset=filtersqp", "QP_solver=BQPD", "hessian_model=LSR1", "quasi_newton_memory_size=20", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) +Optimizer_Uno_filtersqp() = Optimizer(["logger=SILENT", "preset=filtersqp", "QP_solver=BQPD", "hessian_model=LSR1", "quasi_newton_memory_size=20", "max_iterations=10000", "unbounded_objective_threshold=-1e15"]) # This testset runs https://github.com/jump-dev/MINLPTests.jl @testset "MINLPTests" begin diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index 2fbb13248..dea2edcfe 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -373,7 +373,7 @@ void hessian_vector_product(int* int_dimension, const double vector[], const dou } // otherwise, try to perform a Hessian-vector product if possible else if (subproblem->has_hessian_operator()) { - subproblem->compute_hessian_vector_product(subproblem->current_iterate.primals.view(), + subproblem->compute_hessian_vector_product(uno::view(subproblem->current_iterate.primals.data(), dimension), uno::view(vector, dimension), uno::view(result, dimension)); } } \ No newline at end of file From 9dfc9a4cd18875d737f256860ea716a757143cdd Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 19:40:58 +0200 Subject: [PATCH 06/13] Fix in LBFGS --- .../hessian_models/quasi_newton/direct/LBFGSHessian.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp index 78c759d9b..29c301630 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp @@ -90,6 +90,7 @@ namespace uno { if (objective_multiplier != this->fixed_objective_multiplier) { throw std::runtime_error("The quasi-Newton Hessian model was initialized with a different objective multiplier"); } + View model_vector{vector.data(), this->model.number_variables}; // recompute the Hessian representation if the limited memory was updated if (this->hessian_recomputation_required) { @@ -99,7 +100,7 @@ namespace uno { // diagonal contribution δ I for (size_t variable_index: Range(this->model.number_variables)) { - result[variable_index] = this->delta * vector[variable_index]; + result[variable_index] = this->delta * model_vector[variable_index]; } // rank-2 contribution @@ -107,14 +108,14 @@ namespace uno { // work on each column of U (Uᵀ v) for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_U_column = this->U.column(column_index); - const double U_coefficient = -dot(current_U_column, vector); // minus sign for U + const double U_coefficient = -dot(current_U_column, model_vector); // minus sign for U // result += coefficient * current_column blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result.data()); } // work on each column of V (Vᵀ v) for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_V_column = this->V.column(column_index); - const double V_coefficient = dot(current_V_column, vector); // plus sign for V + const double V_coefficient = dot(current_V_column, model_vector); // plus sign for V // result += coefficient * current_column blas1::add(this->model.number_variables, V_coefficient, current_V_column.data(), result.data()); } From cdbe482208a84f47b7abfdc0512ea6045db137ae Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 19:44:19 +0200 Subject: [PATCH 07/13] Fix in L-SR1 --- .../hessian_models/quasi_newton/direct/LBFGSHessian.cpp | 8 ++++---- .../hessian_models/quasi_newton/direct/LSR1Hessian.cpp | 5 +++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp index 29c301630..dd4650d52 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LBFGSHessian.cpp @@ -90,7 +90,7 @@ namespace uno { if (objective_multiplier != this->fixed_objective_multiplier) { throw std::runtime_error("The quasi-Newton Hessian model was initialized with a different objective multiplier"); } - View model_vector{vector.data(), this->model.number_variables}; + View subvector{vector.data(), this->model.number_variables}; // recompute the Hessian representation if the limited memory was updated if (this->hessian_recomputation_required) { @@ -100,7 +100,7 @@ namespace uno { // diagonal contribution δ I for (size_t variable_index: Range(this->model.number_variables)) { - result[variable_index] = this->delta * model_vector[variable_index]; + result[variable_index] = this->delta * subvector[variable_index]; } // rank-2 contribution @@ -108,14 +108,14 @@ namespace uno { // work on each column of U (Uᵀ v) for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_U_column = this->U.column(column_index); - const double U_coefficient = -dot(current_U_column, model_vector); // minus sign for U + const double U_coefficient = -dot(current_U_column, subvector); // minus sign for U // result += coefficient * current_column blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result.data()); } // work on each column of V (Vᵀ v) for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_V_column = this->V.column(column_index); - const double V_coefficient = dot(current_V_column, model_vector); // plus sign for V + const double V_coefficient = dot(current_V_column, subvector); // plus sign for V // result += coefficient * current_column blas1::add(this->model.number_variables, V_coefficient, current_V_column.data(), result.data()); } diff --git a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp index e5d8d5b2b..630f4c469 100644 --- a/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp +++ b/uno/ingredients/hessian_models/quasi_newton/direct/LSR1Hessian.cpp @@ -58,6 +58,7 @@ namespace uno { if (objective_multiplier != this->fixed_objective_multiplier) { throw std::runtime_error("The L-SR1 Hessian model was initialized with a different objective multiplier"); } + View subvector{vector.data(), this->model.number_variables}; // a recomputation of the Hessian representation may be required if (this->hessian_recomputation_required) { @@ -67,7 +68,7 @@ namespace uno { // diagonal contribution δ I for (size_t variable_index: Range(this->model.number_variables)) { - result[variable_index] = this->delta * vector[variable_index]; + result[variable_index] = this->delta * subvector[variable_index]; } // rank-1 contribution: U in R^{n x m} @@ -75,7 +76,7 @@ namespace uno { DEBUG << "U = " << this->U << '\n'; for (size_t column_index: Range(this->number_entries_in_memory)) { const auto current_U_column = this->U.column(column_index); - double U_coefficient = dot(current_U_column, vector) / this->get_correction_column_scaling(column_index); + double U_coefficient = dot(current_U_column, subvector) / this->get_correction_column_scaling(column_index); assert(!std::isnan(U_coefficient)); // result += coefficient * column(P⁻¹) * current_column blas1::add(this->model.number_variables, U_coefficient, current_U_column.data(), result.data()); From d0bcd86333bc3fc0a305d04a2bf64e3085b1b907 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 20:09:56 +0200 Subject: [PATCH 08/13] Cleanup in BQPD --- .../subproblem_solvers/BQPD/BQPDSolver.cpp | 20 +++++++++---------- uno/linear_algebra/View.hpp | 3 +++ 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp index dea2edcfe..ce6ecea30 100644 --- a/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp +++ b/uno/ingredients/subproblem_solvers/BQPD/BQPDSolver.cpp @@ -20,8 +20,8 @@ #define hessian_vector_product FC_GLOBAL(gdotx, GDOTX) extern "C" { - void hessian_vector_product([[maybe_unused]] int *int_dimension, const double vector[], const double ws[], - const int lws[], double result[]); + void hessian_vector_product([[maybe_unused]] int *bqpd_dimension, const double bqpd_vector[], const double ws[], + const int lws[], double bqpd_result[]); // fortran common block used in bqpd/bqpd.f extern struct { @@ -329,13 +329,13 @@ namespace uno { } } // namespace -void hessian_vector_product(int* int_dimension, const double vector[], const double /*ws*/[], const int lws[], double result[]) { - assert(int_dimension != nullptr && *int_dimension >= 0); - const size_t dimension = static_cast(*int_dimension); - - for (size_t i = 0; i < dimension; i++) { - result[i] = 0.; - } +void hessian_vector_product(int* bqpd_dimension, const double bqpd_vector[], const double /*ws*/[], const int lws[], + double bqpd_result[]) { + assert(bqpd_dimension != nullptr && *bqpd_dimension >= 0); + const size_t dimension = static_cast(*bqpd_dimension); + auto result = uno::view(bqpd_result, dimension); + result.fill(0.); + auto vector = uno::view(bqpd_vector, dimension); // retrieve workspace, statistics, and subproblem uno::Statistics* statistics = uno::retrieve_pointer(0, lws); @@ -374,6 +374,6 @@ void hessian_vector_product(int* int_dimension, const double vector[], const dou // otherwise, try to perform a Hessian-vector product if possible else if (subproblem->has_hessian_operator()) { subproblem->compute_hessian_vector_product(uno::view(subproblem->current_iterate.primals.data(), dimension), - uno::view(vector, dimension), uno::view(result, dimension)); + vector, result); } } \ No newline at end of file diff --git a/uno/linear_algebra/View.hpp b/uno/linear_algebra/View.hpp index cf6f01eec..d07627eb7 100644 --- a/uno/linear_algebra/View.hpp +++ b/uno/linear_algebra/View.hpp @@ -4,6 +4,7 @@ #ifndef UNO_VIEW_H #define UNO_VIEW_H +#include #include #include "BLAS.hpp" #include "symbolic/Inverse.hpp" @@ -58,10 +59,12 @@ namespace uno { } [[nodiscard]] T& operator[](size_t index) noexcept { + assert(index < this->size()); return this->pointer[index]; } [[nodiscard]] const T& operator[](size_t index) const noexcept { + assert(index < this->size()); return this->pointer[index]; } From 362d6e4955f4e025c312a974bd4820626e46fc98 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 20:15:02 +0200 Subject: [PATCH 09/13] Removed raw pointers in SymmetricIndefiniteLinearSolver --- uno/ingredients/subproblem_solvers/COOLinearSystem.cpp | 3 +-- uno/ingredients/subproblem_solvers/EQPSolver.cpp | 2 +- uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp | 7 +++---- uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp | 2 +- uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp | 9 ++++----- uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp | 2 +- uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp | 4 ++-- uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp | 2 +- .../SymmetricIndefiniteLinearSolver.hpp | 8 +++----- uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp | 2 +- 10 files changed, 18 insertions(+), 23 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp b/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp index 327c55a8e..207c82bb7 100644 --- a/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp +++ b/uno/ingredients/subproblem_solvers/COOLinearSystem.cpp @@ -36,8 +36,7 @@ namespace uno { this->solution.resize(this->dimension); } - double COOLinearSystem::compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, - const Vector& /*vector*/) const { + double COOLinearSystem::compute_hessian_quadratic_form(const Subproblem& /*subproblem*/, const Vector& /*vector*/) const { return 0.; } } // namespace diff --git a/uno/ingredients/subproblem_solvers/EQPSolver.cpp b/uno/ingredients/subproblem_solvers/EQPSolver.cpp index 7cdeb98bf..17fa18aea 100644 --- a/uno/ingredients/subproblem_solvers/EQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/EQPSolver.cpp @@ -60,7 +60,7 @@ namespace uno { } // solve the linear system - this->linear_solver->solve_indefinite_system(linear_system.solution.data()); + this->linear_solver->solve_indefinite_system(linear_system.solution.view()); if (this->linear_solver->matrix_is_singular()) { direction.status = SubproblemStatus::INFEASIBLE; return; diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp index 4d3b5ed0d..99ed962b0 100644 --- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp +++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.cpp @@ -191,18 +191,17 @@ namespace uno { this->factorization_performed = true; } - void MA27Solver::solve_indefinite_system(double* solution) { + void MA27Solver::solve_indefinite_system(View solution) { assert(this->factorization_performed); int la = static_cast(this->workspace.factor.size()); int liw = static_cast(this->workspace.iw.size()); // copy rhs into solution (overwritten by MA27) - const size_t dimension = static_cast(this->workspace.n); - view(solution, dimension) = this->linear_system.rhs.view(); + solution = this->linear_system.rhs.view(); MA27_linear_solve(&this->workspace.n, this->workspace.factor.data(), &la, this->workspace.iw.data(), &liw, - this->workspace.w.data(), &this->workspace.maxfrt, solution, this->workspace.iw1.data(), &this->workspace.nsteps, + this->workspace.w.data(), &this->workspace.maxfrt, solution.data(), this->workspace.iw1.data(), &this->workspace.nsteps, this->workspace.icntl.data(), this->workspace.info.data()); assert(this->workspace.info[eINFO::IFLAG] == eIFLAG::SUCCESS && "MA27: the linear solve failed"); diff --git a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp index 370638a87..ea7d4d0a5 100644 --- a/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp +++ b/uno/ingredients/subproblem_solvers/MA27/MA27Solver.hpp @@ -44,7 +44,7 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(bool is_matrix_positive_definite) override; - void solve_indefinite_system(double* solution) override; + void solve_indefinite_system(View solution) override; // solve_indefinite_system for multiple RHS is that of the base class using DirectSymmetricIndefiniteLinearSolver::solve_indefinite_system; diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp index f386905b5..960b1aaa4 100644 --- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp +++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.cpp @@ -155,7 +155,7 @@ namespace uno { this->factorization_performed = true; } - void MA57Solver::solve_indefinite_system(double* solution) { + void MA57Solver::solve_indefinite_system(View solution) { assert(this->factorization_performed); // solve @@ -167,17 +167,16 @@ namespace uno { MA57_linear_solve_with_iterative_refinement(&this->workspace.job, &this->workspace.n, &this->workspace.nnz, this->linear_system.matrix_values.data(), this->linear_system.matrix_row_indices.data(), this->linear_system.matrix_column_indices.data(), this->workspace.fact.data(), &this->workspace.lfact, - this->workspace.ifact.data(), &this->workspace.lifact, this->linear_system.rhs.data(), solution, + this->workspace.ifact.data(), &this->workspace.lifact, this->linear_system.rhs.data(), solution.data(), this->workspace.residuals.data(), this->workspace.work.data(), this->workspace.iwork.data(), this->workspace.icntl.data(), this->workspace.cntl.data(), this->workspace.info.data(), this->workspace.rinfo.data()); } else { // copy rhs into solution (overwritten by MA57) - const size_t dimension = static_cast(this->workspace.n); - view(solution, dimension) = this->linear_system.rhs.view(); + solution = this->linear_system.rhs.view(); MA57_linear_solve(&this->workspace.job, &this->workspace.n, this->workspace.fact.data(), &this->workspace.lfact, - this->workspace.ifact.data(), &this->workspace.lifact, &nrhs, solution, &lrhs, + this->workspace.ifact.data(), &this->workspace.lifact, &nrhs, solution.data(), &lrhs, this->workspace.work.data(), &this->workspace.lwork, this->workspace.iwork.data(), this->workspace.icntl.data(), this->workspace.info.data()); } diff --git a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp index 27eb15ee0..5d63e4f63 100644 --- a/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp +++ b/uno/ingredients/subproblem_solvers/MA57/MA57Solver.hpp @@ -50,7 +50,7 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(bool is_matrix_positive_definite) override; - void solve_indefinite_system(double* solution) override; + void solve_indefinite_system(View solution) override; void solve_indefinite_system(const double* rhs, double* solution, size_t number_of_rhs) override; [[nodiscard]] Inertia get_inertia() const override; diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp index 38d974cde..0f65d7df5 100644 --- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp +++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.cpp @@ -100,8 +100,8 @@ namespace uno { this->factorization_performed = true; } - void MUMPSSolver::solve_indefinite_system(double* solution) { - return this->solve_indefinite_system(this->linear_system.rhs.data(), solution, 1); + void MUMPSSolver::solve_indefinite_system(View solution) { + return this->solve_indefinite_system(this->linear_system.rhs.data(), solution.data(), 1); } void MUMPSSolver::solve_indefinite_system(const double* rhs, double* solution, size_t number_of_rhs) { diff --git a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp index 24035daa0..8eaf31e9a 100644 --- a/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp +++ b/uno/ingredients/subproblem_solvers/MUMPS/MUMPSSolver.hpp @@ -31,7 +31,7 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(bool is_matrix_positive_definite) override; - void solve_indefinite_system(double* solution) override; + void solve_indefinite_system(View solution) override; void solve_indefinite_system(const double* rhs, double* solution, size_t number_of_rhs) override; [[nodiscard]] Inertia get_inertia() const override; diff --git a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp index 55d292dfb..a1ddb4d9e 100644 --- a/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp +++ b/uno/ingredients/subproblem_solvers/SymmetricIndefiniteLinearSolver.hpp @@ -6,11 +6,9 @@ #include #include "LinearSystem.hpp" +#include "linear_algebra/View.hpp" namespace uno { - // forward declarations - class LinearSystem; - template class SymmetricIndefiniteLinearSolver { public: @@ -21,7 +19,7 @@ namespace uno { // solve with a single right-hand side: the right-hand side is read from the linear system // (get_linear_system().rhs) and the solution is written to `solution` - virtual void solve_indefinite_system(ElementType* solution) = 0; + virtual void solve_indefinite_system(View solution) = 0; // solve with `number_of_rhs` right-hand sides stored column-major (each column has length // get_linear_system().dimension): `rhs` is the block of right-hand sides and `solution` receives @@ -36,7 +34,7 @@ namespace uno { const auto rhs_column = view(rhs, column_offset, column_offset + linear_system.dimension); // copy the column into linear_system.rhs, then solve into the matching solution column linear_system.rhs.view() = rhs_column; - double* rhs_solution = solution + column_offset; + auto rhs_solution = view(solution, column_offset, column_offset + linear_system.dimension); this->solve_indefinite_system(rhs_solution); } } diff --git a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp index b76fd05f1..fdc746bfa 100644 --- a/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp +++ b/uno/ingredients/subproblem_solvers/WoodburyEQPSolver.cpp @@ -63,7 +63,7 @@ namespace uno { // solve the linear system with only the diagonal part and store the result in solution_diagonal_part Vector solution_diagonal_part(subproblem.number_variables + subproblem.number_constraints); - this->linear_solver->solve_indefinite_system(solution_diagonal_part.data()); + this->linear_solver->solve_indefinite_system(solution_diagonal_part.view()); if (this->linear_solver->matrix_is_singular()) { direction.status = SubproblemStatus::INFEASIBLE; return; From c6471c76344114ab09de85d0c98900657ec4ae86 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 20:20:30 +0200 Subject: [PATCH 10/13] Fixes --- uno/linear_algebra/BLAS.hpp | 1 + unotest/functional_tests/MA27SolverTests.cpp | 2 +- unotest/functional_tests/MA57SolverTests.cpp | 2 +- unotest/functional_tests/MUMPSSolverTests.cpp | 2 +- unotest/functional_tests/SSIDSSolverTests.cpp | 2 +- 5 files changed, 5 insertions(+), 4 deletions(-) diff --git a/uno/linear_algebra/BLAS.hpp b/uno/linear_algebra/BLAS.hpp index 8a3b392f1..beb0221d8 100644 --- a/uno/linear_algebra/BLAS.hpp +++ b/uno/linear_algebra/BLAS.hpp @@ -4,6 +4,7 @@ #ifndef UNO_BLAS_H #define UNO_BLAS_H +#include #include #include "fortran_interface.h" #define dcopy FC_GLOBAL_(dcopy, DCOPY) diff --git a/unotest/functional_tests/MA27SolverTests.cpp b/unotest/functional_tests/MA27SolverTests.cpp index a77bc5140..3311aaa67 100644 --- a/unotest/functional_tests/MA27SolverTests.cpp +++ b/unotest/functional_tests/MA27SolverTests.cpp @@ -24,7 +24,7 @@ TEST(MA27Solver, SystemSize5) { solver.do_numerical_factorization(false); Vector result(dimension); result.fill(0.); - solver.solve_indefinite_system(result.data()); + solver.solve_indefinite_system(result.view()); const std::array reference{1., 2., 3., 4., 5.}; const double tolerance = 1e-8; diff --git a/unotest/functional_tests/MA57SolverTests.cpp b/unotest/functional_tests/MA57SolverTests.cpp index 8f7ed21eb..e35517dfa 100644 --- a/unotest/functional_tests/MA57SolverTests.cpp +++ b/unotest/functional_tests/MA57SolverTests.cpp @@ -24,7 +24,7 @@ TEST(MA57Solver, SystemSize5) { solver.do_numerical_factorization(false); Vector result(dimension); result.fill(0.); - solver.solve_indefinite_system(result.data()); + solver.solve_indefinite_system(result.view()); const std::array reference{1., 2., 3., 4., 5.}; const double tolerance = 1e-8; diff --git a/unotest/functional_tests/MUMPSSolverTests.cpp b/unotest/functional_tests/MUMPSSolverTests.cpp index 4d29b0e72..75c98ee8c 100644 --- a/unotest/functional_tests/MUMPSSolverTests.cpp +++ b/unotest/functional_tests/MUMPSSolverTests.cpp @@ -24,7 +24,7 @@ TEST(MUMPSSolver, SystemSize5) { solver.do_numerical_factorization(false); Vector result(dimension); result.fill(0.); - solver.solve_indefinite_system(result.data()); + solver.solve_indefinite_system(result.view()); const std::array reference{1., 2., 3., 4., 5.}; const double tolerance = 1e-8; diff --git a/unotest/functional_tests/SSIDSSolverTests.cpp b/unotest/functional_tests/SSIDSSolverTests.cpp index 2119d93b5..9234d82ae 100644 --- a/unotest/functional_tests/SSIDSSolverTests.cpp +++ b/unotest/functional_tests/SSIDSSolverTests.cpp @@ -24,7 +24,7 @@ TEST(SSIDSSolver, SystemSize5) { solver.do_numerical_factorization(false); Vector result(dimension); result.fill(0.); - solver.solve_indefinite_system(result.data()); + solver.solve_indefinite_system(result.view()); const std::array reference{1., 2., 3., 4., 5.}; const double tolerance = 1e-8; From f736c6e871e0d1657ac1c5543f246d7a540bee5f Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Thu, 11 Jun 2026 20:26:58 +0200 Subject: [PATCH 11/13] Fixed SSIDS --- uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.cpp | 7 +++---- uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.hpp | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.cpp b/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.cpp index 9243eb2ad..0cab71888 100644 --- a/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.cpp +++ b/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.cpp @@ -42,14 +42,13 @@ namespace uno { this->factorization_performed = true; } - void SSIDSSolver::solve_indefinite_system(double* solution) { + void SSIDSSolver::solve_indefinite_system(View solution) { assert(this->factorization_performed); // copy rhs into solution (overwritten by SSIDS) - const size_t dimension = static_cast(this->workspace.n); - view(solution, dimension) = this->linear_system.rhs.view(); + solution = this->linear_system.rhs.view(); - spral_ssids_solve1(0, solution, this->workspace.akeep, this->workspace.fkeep, &this->workspace.options, + spral_ssids_solve1(0, solution.data(), this->workspace.akeep, this->workspace.fkeep, &this->workspace.options, &this->workspace.inform); if (this->workspace.inform.flag < 0) { spral_ssids_free(&this->workspace.akeep, &this->workspace.fkeep); diff --git a/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.hpp b/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.hpp index ff7aa5619..a1c0f17d8 100644 --- a/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.hpp +++ b/uno/ingredients/subproblem_solvers/SSIDS/SSIDSSolver.hpp @@ -28,7 +28,7 @@ namespace uno { void do_symbolic_analysis() override; void do_numerical_factorization(bool is_matrix_positive_definite) override; - void solve_indefinite_system(double* solution) override; + void solve_indefinite_system(View solution) override; void solve_indefinite_system(const double* rhs, double* solution, size_t number_of_rhs) override; [[nodiscard]] Inertia get_inertia() const override; From 2643250c578357d13e631a4b60c1b4c9e58a2ba5 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 12 Jun 2026 14:54:27 +0200 Subject: [PATCH 12/13] Progress --- interfaces/AMPL/AMPLModel.cpp | 5 +++-- interfaces/AMPL/AMPLModel.hpp | 4 ++-- interfaces/C/Uno_C_API.cpp | 8 ++++---- interfaces/Python/cpp_classes/PythonModel.cpp | 16 ++++++++-------- interfaces/Python/cpp_classes/PythonModel.hpp | 4 ++-- .../FeasibilityRestoration.cpp | 2 +- .../relaxed_problems/l1RelaxedProblem.cpp | 9 +++++---- .../relaxed_problems/l1RelaxedProblem.hpp | 5 +++-- .../PrimalDualInteriorPointProblem.cpp | 4 ++-- .../PrimalDualInteriorPointProblem.hpp | 5 +++-- uno/ingredients/subproblem/Subproblem.cpp | 6 +++--- uno/model/BoundRelaxedModel.hpp | 4 ++-- uno/model/FixedBoundsConstraintsModel.cpp | 7 ++++--- uno/model/FixedBoundsConstraintsModel.hpp | 4 ++-- .../HomogeneousEqualityConstrainedModel.cpp | 7 ++++--- .../HomogeneousEqualityConstrainedModel.hpp | 5 +++-- uno/model/Model.cpp | 2 +- uno/model/Model.hpp | 4 ++-- uno/optimization/Evaluations.cpp | 6 ++++-- uno/optimization/Evaluations.hpp | 4 ++-- uno/optimization/OptimizationProblem.cpp | 4 ++-- uno/optimization/OptimizationProblem.hpp | 4 ++-- 22 files changed, 64 insertions(+), 55 deletions(-) diff --git a/interfaces/AMPL/AMPLModel.cpp b/interfaces/AMPL/AMPLModel.cpp index 2767311f1..59456111a 100644 --- a/interfaces/AMPL/AMPLModel.cpp +++ b/interfaces/AMPL/AMPLModel.cpp @@ -217,11 +217,12 @@ namespace uno { ++this->number_model_evaluations.hessian; } - void AMPLModel::compute_jacobian_vector_product(const double* /*x*/, const double* /*vector*/, double* /*result*/) const { + void AMPLModel::compute_jacobian_vector_product(View /*x*/, View /*vector*/, View /*result*/) const { throw std::runtime_error("AMPLModel::compute_jacobian_vector_product not implemented"); } - void AMPLModel::compute_jacobian_transposed_vector_product(const double* /*x*/, const double* /*vector*/, double* /*result*/) const { + void AMPLModel::compute_jacobian_transposed_vector_product(View /*x*/, View /*vector*/, + View /*result*/) const { throw std::runtime_error("AMPLModel::compute_jacobian_transposed_vector_product not implemented"); } diff --git a/interfaces/AMPL/AMPLModel.hpp b/interfaces/AMPL/AMPLModel.hpp index 331b3b591..2b289acaf 100644 --- a/interfaces/AMPL/AMPLModel.hpp +++ b/interfaces/AMPL/AMPLModel.hpp @@ -55,8 +55,8 @@ namespace uno { double* hessian_values) const override; // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; + void compute_jacobian_vector_product(View x, View vector, View result) const override; + void compute_jacobian_transposed_vector_product(View x, View vector, View result) const override; void compute_hessian_vector_product(View x, View vector, double objective_multiplier, const Vector& multipliers, View result) const override; diff --git a/interfaces/C/Uno_C_API.cpp b/interfaces/C/Uno_C_API.cpp index c08df82a3..6de4672c3 100644 --- a/interfaces/C/Uno_C_API.cpp +++ b/interfaces/C/Uno_C_API.cpp @@ -184,10 +184,10 @@ class UnoModel: public Model { } } - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override { + void compute_jacobian_vector_product(View x, View vector, View result) const override { if (this->user_model.jacobian_operator != nullptr) { const uno_int return_code = this->user_model.jacobian_operator(this->user_model.number_variables, - this->user_model.number_constraints, x, true, vector, result, this->user_model.user_data); + this->user_model.number_constraints, x.data(), true, vector.data(), result.data(), this->user_model.user_data); if (0 < return_code) { throw GradientEvaluationError(); } @@ -197,10 +197,10 @@ class UnoModel: public Model { } } - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override { + void compute_jacobian_transposed_vector_product(View x, View vector, View result) const override { if (this->user_model.jacobian_transposed_operator != nullptr) { const uno_int return_code = this->user_model.jacobian_transposed_operator(this->user_model.number_variables, - this->user_model.number_constraints, x, true, vector, result, this->user_model.user_data); + this->user_model.number_constraints, x.data(), true, vector.data(), result.data(), this->user_model.user_data); if (0 < return_code) { throw GradientEvaluationError(); } diff --git a/interfaces/Python/cpp_classes/PythonModel.cpp b/interfaces/Python/cpp_classes/PythonModel.cpp index c6d5ac6a4..f157dc97a 100644 --- a/interfaces/Python/cpp_classes/PythonModel.cpp +++ b/interfaces/Python/cpp_classes/PythonModel.cpp @@ -181,11 +181,11 @@ namespace uno { } } - void PythonModel::compute_jacobian_vector_product(const double* x, const double* vector, double* result) const { + void PythonModel::compute_jacobian_vector_product(View x, View vector, View result) const { if (this->user_model.jacobian_operator.has_value()) { - const auto x_py = to_const_array(x, this->number_variables); - const auto vector_py = to_const_array(vector, this->number_variables); - auto result_py = to_array(result, this->number_constraints); + const auto x_py = to_const_array(x.data(), this->number_variables); + const auto vector_py = to_const_array(vector.data(), this->number_variables); + auto result_py = to_array(result.data(), this->number_constraints); // evaluate Jacobian-vector product try { @@ -200,11 +200,11 @@ namespace uno { } } - void PythonModel::compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const { + void PythonModel::compute_jacobian_transposed_vector_product(View x, View vector, View result) const { if (this->user_model.jacobian_transposed_operator.has_value()) { - const auto x_py = to_const_array(x, this->number_variables); - const auto vector_py = to_const_array(vector, this->number_constraints); - auto result_py = to_array(result, this->number_variables); + const auto x_py = to_const_array(x.data(), this->number_variables); + const auto vector_py = to_const_array(vector.data(), this->number_constraints); + auto result_py = to_array(result.data(), this->number_variables); // evaluate Jacobian^T-vector product try { diff --git a/interfaces/Python/cpp_classes/PythonModel.hpp b/interfaces/Python/cpp_classes/PythonModel.hpp index bf1083609..b3e05502e 100644 --- a/interfaces/Python/cpp_classes/PythonModel.hpp +++ b/interfaces/Python/cpp_classes/PythonModel.hpp @@ -45,8 +45,8 @@ namespace uno { double* hessian_values) const override; // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; + void compute_jacobian_vector_product(View x, View vector, View result) const override; + void compute_jacobian_transposed_vector_product(View x, View vector, View result) const override; void compute_hessian_vector_product(View x, View vector, double objective_multiplier, const Vector& multipliers, View result) const override; diff --git a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp index 7b9dfc8fb..651a8857d 100644 --- a/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/FeasibilityRestoration.cpp @@ -194,7 +194,7 @@ namespace uno { // compute the linearized constraint violation // TODO preallocate Vector result(model.number_constraints); - current_evaluations.compute_jacobian_vector_product(model, direction.primals.data(), result.data()); + current_evaluations.compute_jacobian_vector_product(model, direction.primals.view(), result.view()); const double trial_linearized_constraint_violation = model.constraint_violation(current_evaluations.constraints + step_length * result, this->residual_norm); return (trial_linearized_constraint_violation <= this->linear_feasibility_tolerance); diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp index 133eab763..eb4f18d79 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.cpp @@ -178,8 +178,8 @@ namespace uno { // ∇c(x_k) λ_k evaluations.evaluate_jacobian(this->model, iterate.primals); - evaluations.compute_jacobian_transposed_vector_product(this->model, iterate.multipliers.constraints.data(), - lagrangian_gradient.data()); + evaluations.compute_jacobian_transposed_vector_product(this->model, iterate.multipliers.constraints.view(), + lagrangian_gradient.view()); lagrangian_gradient.scale(-1.); // z_k @@ -238,7 +238,8 @@ namespace uno { } } - void l1RelaxedProblem::compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const { + void l1RelaxedProblem::compute_jacobian_vector_product(View vector, View result, + const Evaluations& evaluations) const { evaluations.compute_jacobian_vector_product(this->model, vector, result); // add the contribution of the elastic variables @@ -250,7 +251,7 @@ namespace uno { } } - void l1RelaxedProblem::compute_jacobian_transposed_vector_product(const double* vector, double* result, + void l1RelaxedProblem::compute_jacobian_transposed_vector_product(View vector, View result, const Evaluations& evaluations) const { evaluations.compute_jacobian_transposed_vector_product(this->model, vector, result); diff --git a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp index b39a6e742..0877b7483 100644 --- a/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp +++ b/uno/ingredients/constraint_relaxation_strategies/relaxed_problems/l1RelaxedProblem.hpp @@ -44,8 +44,9 @@ namespace uno { const Multipliers& multipliers, double* hessian_values) const override; // linear operators - void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; - void compute_jacobian_transposed_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; + void compute_jacobian_vector_product(View vector, View result, const Evaluations& evaluations) const override; + void compute_jacobian_transposed_vector_product(View vector, View result, + const Evaluations& evaluations) const override; void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, const Multipliers& multipliers, View result) const override; diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp index 2b4409506..438691a8a 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.cpp @@ -232,12 +232,12 @@ namespace uno { } } - void PrimalDualInteriorPointProblem::compute_jacobian_vector_product(const double* vector, double* result, + void PrimalDualInteriorPointProblem::compute_jacobian_vector_product(View vector, View result, const Evaluations& evaluations) const { this->inner.compute_jacobian_vector_product(vector, result, evaluations); } - void PrimalDualInteriorPointProblem::compute_jacobian_transposed_vector_product(const double* vector, double* result, + void PrimalDualInteriorPointProblem::compute_jacobian_transposed_vector_product(View vector, View result, const Evaluations& evaluations) const { this->inner.compute_jacobian_transposed_vector_product(vector, result, evaluations); } diff --git a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp index 42e0dc905..f74af0dd6 100644 --- a/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp +++ b/uno/ingredients/inequality_handling_methods/interior_point_methods/barrier_problems/PrimalDualInteriorPointProblem.hpp @@ -44,8 +44,9 @@ namespace uno { const Multipliers& multipliers, double* hessian_values) const override; // linear operators - void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; - void compute_jacobian_transposed_vector_product(const double* vector, double* result, const Evaluations& evaluations) const override; + void compute_jacobian_vector_product(View vector, View result, const Evaluations& evaluations) const override; + void compute_jacobian_transposed_vector_product(View vector, View result, + const Evaluations& evaluations) const override; void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, const Multipliers& multipliers, View result) const override; diff --git a/uno/ingredients/subproblem/Subproblem.cpp b/uno/ingredients/subproblem/Subproblem.cpp index efd128768..0443be3a1 100644 --- a/uno/ingredients/subproblem/Subproblem.cpp +++ b/uno/ingredients/subproblem/Subproblem.cpp @@ -135,8 +135,8 @@ namespace uno { rhs.fill(0.); // -Jacobian^T-multipliers product - this->problem.compute_jacobian_transposed_vector_product(this->current_iterate.multipliers.constraints.data(), - rhs.data(), evaluations); + this->problem.compute_jacobian_transposed_vector_product(this->current_iterate.multipliers.constraints.view(), + rhs.view(), evaluations); rhs.scale(-1.); // objective gradient @@ -271,7 +271,7 @@ namespace uno { const double current_constraint_violation = model.constraint_violation(current_evaluations.constraints, norm); // TODO preallocate Vector result(model.number_constraints); - current_evaluations.compute_jacobian_vector_product(model, primal_direction.data(), result.data()); + current_evaluations.compute_jacobian_vector_product(model, primal_direction.view(), result.view()); const double trial_linearized_constraint_violation = model.constraint_violation(current_evaluations.constraints + step_length * result, norm); return current_constraint_violation - trial_linearized_constraint_violation; diff --git a/uno/model/BoundRelaxedModel.hpp b/uno/model/BoundRelaxedModel.hpp index 8ce46fee7..15f6636f2 100644 --- a/uno/model/BoundRelaxedModel.hpp +++ b/uno/model/BoundRelaxedModel.hpp @@ -69,11 +69,11 @@ namespace uno { this->model.evaluate_lagrangian_hessian(x, objective_multiplier, multipliers, hessian_values); } - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override { + void compute_jacobian_vector_product(View x, View vector, View result) const override { this->model.compute_jacobian_vector_product(x, vector, result); } - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override { + void compute_jacobian_transposed_vector_product(View x, View vector, View result) const override { this->model.compute_jacobian_transposed_vector_product(x, vector, result); } diff --git a/uno/model/FixedBoundsConstraintsModel.cpp b/uno/model/FixedBoundsConstraintsModel.cpp index 8e407776b..cf625b20f 100644 --- a/uno/model/FixedBoundsConstraintsModel.cpp +++ b/uno/model/FixedBoundsConstraintsModel.cpp @@ -110,7 +110,8 @@ namespace uno { } // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void FixedBoundsConstraintsModel::compute_jacobian_vector_product(const double* x, const double* vector, double* result) const { + void FixedBoundsConstraintsModel::compute_jacobian_vector_product(View x, View vector, + View result) const { this->model.compute_jacobian_vector_product(x, vector, result); // add the contributions of the fixed variables @@ -121,8 +122,8 @@ namespace uno { } } - void FixedBoundsConstraintsModel::compute_jacobian_transposed_vector_product(const double* x, const double* vector, - double* result) const { + void FixedBoundsConstraintsModel::compute_jacobian_transposed_vector_product(View x, View vector, + View result) const { this->model.compute_jacobian_transposed_vector_product(x, vector, result); // add the contributions of the fixed variables diff --git a/uno/model/FixedBoundsConstraintsModel.hpp b/uno/model/FixedBoundsConstraintsModel.hpp index acb136c13..59493b39c 100644 --- a/uno/model/FixedBoundsConstraintsModel.hpp +++ b/uno/model/FixedBoundsConstraintsModel.hpp @@ -46,8 +46,8 @@ namespace uno { double* hessian_values) const override; // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; + void compute_jacobian_vector_product(View x, View vector, View result) const override; + void compute_jacobian_transposed_vector_product(View x, View vector, View result) const override; void compute_hessian_vector_product(View x, View vector, double objective_multiplier, const Vector& multipliers, View result) const override; diff --git a/uno/model/HomogeneousEqualityConstrainedModel.cpp b/uno/model/HomogeneousEqualityConstrainedModel.cpp index ecefd7ccf..e46079a2c 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.cpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.cpp @@ -115,7 +115,8 @@ namespace uno { } // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void HomogeneousEqualityConstrainedModel::compute_jacobian_vector_product(const double* x, const double* vector, double* result) const { + void HomogeneousEqualityConstrainedModel::compute_jacobian_vector_product(View x, View vector, + View result) const { this->model.compute_jacobian_vector_product(x, vector, result); // add the slack contributions @@ -124,8 +125,8 @@ namespace uno { } } - void HomogeneousEqualityConstrainedModel::compute_jacobian_transposed_vector_product(const double* x, const double* vector, - double* result) const { + void HomogeneousEqualityConstrainedModel::compute_jacobian_transposed_vector_product(View x, View vector, + View result) const { this->model.compute_jacobian_transposed_vector_product(x, vector, result); // add the slack contributions diff --git a/uno/model/HomogeneousEqualityConstrainedModel.hpp b/uno/model/HomogeneousEqualityConstrainedModel.hpp index 2cad620f8..83b038118 100644 --- a/uno/model/HomogeneousEqualityConstrainedModel.hpp +++ b/uno/model/HomogeneousEqualityConstrainedModel.hpp @@ -43,8 +43,9 @@ namespace uno { double* hessian_values) const override; // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products - void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const override; - void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const override; + void compute_jacobian_vector_product(View x, View vector, View result) const override; + void compute_jacobian_transposed_vector_product(View x, View vector, + View result) const override; void compute_hessian_vector_product(View x, View vector, double objective_multiplier, const Vector& multipliers, View result) const override; diff --git a/uno/model/Model.cpp b/uno/model/Model.cpp index e2efaa739..6250f00ec 100644 --- a/uno/model/Model.cpp +++ b/uno/model/Model.cpp @@ -52,7 +52,7 @@ namespace uno { if (0 < this->number_constraints) { // TODO test whether λ_k != 0 evaluations.evaluate_jacobian(*this, primals); - evaluations.compute_jacobian_transposed_vector_product(*this, multipliers.constraints.data(), lagrangian_gradient.data()); + evaluations.compute_jacobian_transposed_vector_product(*this, multipliers.constraints.view(), lagrangian_gradient.view()); lagrangian_gradient.scale(-1.); } diff --git a/uno/model/Model.hpp b/uno/model/Model.hpp index cee1f2c58..390c3a516 100644 --- a/uno/model/Model.hpp +++ b/uno/model/Model.hpp @@ -73,8 +73,8 @@ namespace uno { // linear operators for Jacobian-, Jacobian^T-, and Hessian-vector products // here we use pointers, since the vector and the result may be provided by a low-level subproblem solver - virtual void compute_jacobian_vector_product(const double* x, const double* vector, double* result) const = 0; - virtual void compute_jacobian_transposed_vector_product(const double* x, const double* vector, double* result) const = 0; + virtual void compute_jacobian_vector_product(View x, View vector, View result) const = 0; + virtual void compute_jacobian_transposed_vector_product(View x, View vector, View result) const = 0; virtual void compute_hessian_vector_product(View x, View vector, double objective_multiplier, const Vector& multipliers, View result) const = 0; diff --git a/uno/optimization/Evaluations.cpp b/uno/optimization/Evaluations.cpp index 6ea21567f..58ba2d395 100644 --- a/uno/optimization/Evaluations.cpp +++ b/uno/optimization/Evaluations.cpp @@ -67,7 +67,8 @@ namespace uno { } // y = J v, where J has dimensions (m, n), v has dimensions (n, 1), and y has dimensions (m, 1) - void Evaluations::compute_jacobian_vector_product([[maybe_unused]] const Model& model, const double* vector, double* result) const { + void Evaluations::compute_jacobian_vector_product([[maybe_unused]] const Model& model, View vector, + View result) const { const size_t number_jacobian_nonzeros = this->jacobian_sparsity->row_indices.size(); for (size_t nonzero_index: Range(number_jacobian_nonzeros)) { const size_t constraint_index = static_cast(this->jacobian_sparsity->row_indices[nonzero_index]); @@ -82,7 +83,8 @@ namespace uno { } // x = Jᵀv, where J has dimensions (m, n), v has dimensions (m, 1), and x has dimensions (n, 1) - void Evaluations::compute_jacobian_transposed_vector_product([[maybe_unused]] const Model& model, const double* vector, double* result) const { + void Evaluations::compute_jacobian_transposed_vector_product([[maybe_unused]] const Model& model, View vector, + View result) const { const size_t number_jacobian_nonzeros = this->jacobian_sparsity->row_indices.size(); for (size_t nonzero_index: Range(number_jacobian_nonzeros)) { const size_t constraint_index = static_cast(this->jacobian_sparsity->row_indices[nonzero_index]); diff --git a/uno/optimization/Evaluations.hpp b/uno/optimization/Evaluations.hpp index 91e322aaf..3d29704ab 100644 --- a/uno/optimization/Evaluations.hpp +++ b/uno/optimization/Evaluations.hpp @@ -37,8 +37,8 @@ namespace uno { void evaluate_objective_gradient(const Model& model, const Vector& primals); void evaluate_jacobian(const Model& model, const Vector& primals); - void compute_jacobian_vector_product(const Model& model, const double* vector, double* result) const; - void compute_jacobian_transposed_vector_product(const Model& model, const double* vector, double* result) const; + void compute_jacobian_vector_product(const Model& model, View vector, View result) const; + void compute_jacobian_transposed_vector_product(const Model& model, View vector, View result) const; // reset the evaluation flags void reset(); diff --git a/uno/optimization/OptimizationProblem.cpp b/uno/optimization/OptimizationProblem.cpp index 1b18c663b..9d58eefcc 100644 --- a/uno/optimization/OptimizationProblem.cpp +++ b/uno/optimization/OptimizationProblem.cpp @@ -104,12 +104,12 @@ namespace uno { multipliers.constraints, hessian_values); } - void OptimizationProblem::compute_jacobian_vector_product(const double* vector, double* result, + void OptimizationProblem::compute_jacobian_vector_product(View vector, View result, const Evaluations& evaluations) const { evaluations.compute_jacobian_vector_product(this->model, vector, result); } - void OptimizationProblem::compute_jacobian_transposed_vector_product(const double* vector, double* result, + void OptimizationProblem::compute_jacobian_transposed_vector_product(View vector, View result, const Evaluations& evaluations) const { evaluations.compute_jacobian_transposed_vector_product(this->model, vector, result); } diff --git a/uno/optimization/OptimizationProblem.hpp b/uno/optimization/OptimizationProblem.hpp index 0c2e1c85c..3a330c165 100644 --- a/uno/optimization/OptimizationProblem.hpp +++ b/uno/optimization/OptimizationProblem.hpp @@ -68,8 +68,8 @@ namespace uno { const Vector& primal_variables, const Multipliers& multipliers, double* hessian_values) const; // linear operators - virtual void compute_jacobian_vector_product(const double* vector, double* result, const Evaluations& evaluations) const; - virtual void compute_jacobian_transposed_vector_product(const double* vector, double* result, + virtual void compute_jacobian_vector_product(View vector, View result, const Evaluations& evaluations) const; + virtual void compute_jacobian_transposed_vector_product(View vector, View result, const Evaluations& evaluations) const; virtual void compute_hessian_vector_product(HessianModel& hessian_model, View x, View vector, const Multipliers& multipliers, View result) const; From d9bb038a2fdb60fb54def48d5ac8fd2d32861332 Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Fri, 19 Jun 2026 12:08:40 +0200 Subject: [PATCH 13/13] Check the sizes --- uno/linear_algebra/View.hpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/uno/linear_algebra/View.hpp b/uno/linear_algebra/View.hpp index d07627eb7..57ed8e5cb 100644 --- a/uno/linear_algebra/View.hpp +++ b/uno/linear_algebra/View.hpp @@ -6,6 +6,7 @@ #include #include +#include #include "BLAS.hpp" #include "symbolic/Inverse.hpp" #include "symbolic/Multiplication.hpp" @@ -37,6 +38,9 @@ namespace uno { View(View&& other) = default; View& operator=(const View& other) { if (&other != this) { + if (this->size() != other.size()) { + throw std::runtime_error("The two views have different sizes"); + } blas1::copy(this->size(), other.data(), this->data()); } return *this;