Model: removed flag fixed_hessian_sparsity + cleanup

bindings/AMPL/AMPLModel.cpp

@@ -79,8 +79,8 @@ namespace uno {
       this->linear_constraints.reserve(this->number_constraints);
       this->generate_constraints();
 
-      // compute number of nonzeros in the Lagrangian Hessian
-      this->set_number_hessian_nonzeros();
+      // compute sparsity pattern and number of nonzeros of Lagrangian Hessian
+      this->compute_lagrangian_hessian_sparsity();
    }
 
    AMPLModel::~AMPLModel() {
@@ -199,13 +199,18 @@ namespace uno {
       }
    }
 
-   void AMPLModel::set_number_hessian_nonzeros() {
+   void AMPLModel::compute_lagrangian_hessian_sparsity() {
       // compute the maximum number of nonzero elements, provided that all multipliers are non-zero
       // int (*Sphset) (ASL*, SputInfo**, int nobj, int ow, int y, int uptri);
       const int objective_number = -1;
       const int upper_triangular = 1;
       this->number_asl_hessian_nonzeros = static_cast<size_t>((*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, 1, 1, upper_triangular));
       this->asl_hessian.reserve(this->number_asl_hessian_nonzeros);
+
+      // sparsity pattern
+      [[maybe_unused]] const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;
+      // check that the column pointers are sorted in increasing order
+      assert(in_increasing_order(asl_column_start, this->number_variables + 1) && "AMPLModel::evaluate_lagrangian_hessian: column starts are not ordered");
    }
 
    size_t AMPLModel::number_objective_gradient_nonzeros() const {
@@ -252,49 +257,20 @@ namespace uno {
       return this->upper_bounded_variables_collection;
    }
 
-   bool are_all_zeros(const Vector<double>& multipliers) {
-      return std::all_of(multipliers.begin(), multipliers.end(), [](double xj) {
-         return xj == 0.;
-      });
-   }
-
-   size_t AMPLModel::compute_hessian_number_nonzeros(double objective_multiplier, const Vector<double>& multipliers) const {
-      // compute the sparsity
-      const int objective_number = -1;
-      const int upper_triangular = 1;
-      const bool all_zeros_multipliers = are_all_zeros(multipliers);
-      int number_nonzeros = (*(this->asl)->p.Sphset)(this->asl, nullptr, objective_number, (objective_multiplier != 0.),
-            not all_zeros_multipliers, upper_triangular);
-      return static_cast<size_t>(number_nonzeros);
-   }
-
    void AMPLModel::evaluate_lagrangian_hessian(const Vector<double>& x, double objective_multiplier, const Vector<double>& multipliers,
          SymmetricMatrix<size_t, double>& hessian) const {
+      assert(hessian.capacity() >= this->number_asl_hessian_nonzeros);
+
       // register the vector of variables
       (*(this->asl)->p.Xknown)(this->asl, const_cast<double*>(x.data()), nullptr);
 
-      // scale by the objective sign
-      objective_multiplier *= this->objective_sign;
-
-      // compute the number of nonzeros
-      [[maybe_unused]] const size_t number_nonzeros = this->fixed_hessian_sparsity ? this->number_asl_hessian_nonzeros :
-                                                      this->compute_hessian_number_nonzeros(objective_multiplier, multipliers);
-      assert(hessian.capacity() >= number_nonzeros);
-
       // evaluate the Hessian: store the matrix in a preallocated array this->asl_hessian
       const int objective_number = -1;
+      objective_multiplier *= this->objective_sign;
       // flip the signs of the multipliers: in AMPL, the Lagrangian is f + lambda.g, while Uno uses f - lambda.g
       this->multipliers_with_flipped_sign = -multipliers;
-      if (this->fixed_hessian_sparsity) {
-         (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, &objective_multiplier,
-               const_cast<double*>(this->multipliers_with_flipped_sign.data()));
-      }
-      else {
-         double* objective_multiplier_pointer = (objective_multiplier != 0.) ? &objective_multiplier : nullptr;
-         const bool all_zeros_multipliers = are_all_zeros(multipliers);
-         (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, objective_multiplier_pointer,
-               all_zeros_multipliers ? nullptr : const_cast<double*>(this->multipliers_with_flipped_sign.data()));
-      }
+      (*(this->asl)->p.Sphes)(this->asl, nullptr, const_cast<double*>(this->asl_hessian.data()), objective_number, &objective_multiplier,
+            const_cast<double*>(this->multipliers_with_flipped_sign.data()));
 
       // generate the sparsity pattern in the right sparse format
       const fint* asl_column_start = this->asl->i.sputinfo_->hcolstarts;

bindings/AMPL/AMPLModel.hpp

@@ -107,8 +107,7 @@ namespace uno {
       void generate_variables();
       void generate_constraints();
 
-      void set_number_hessian_nonzeros();
-      [[nodiscard]] size_t compute_hessian_number_nonzeros(double objective_multiplier, const Vector<double>& multipliers) const;
+      void compute_lagrangian_hessian_sparsity();
       static void determine_bounds_types(const std::vector<double>& lower_bounds, const std::vector<double>& upper_bounds, std::vector<BoundType>& status);
    };
 

uno/linear_algebra/SymmetricIndefiniteLinearSystem.hpp

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

uno/model/Model.hpp

@@ -44,9 +44,6 @@ namespace uno {
       const size_t number_constraints; /*!< Number of constraints */
       const double objective_sign; /*!< Sign of the objective function (1: minimization, -1: maximization) */
 
-      // Hessian
-      const bool fixed_hessian_sparsity{true};
-
       [[nodiscard]] virtual double evaluate_objective(const Vector<double>& x) const = 0;
       virtual void evaluate_objective_gradient(const Vector<double>& x, SparseVector<double>& gradient) const = 0;
       virtual void evaluate_constraints(const Vector<double>& x, std::vector<double>& constraints) const = 0;