Moved Hessian-vector product gdotx from wdotd.f (Fortran) to BQPDSolv…

…er.cpp (C++)

CMakeLists.txt

@@ -142,7 +142,7 @@ find_library(BQPD bqpd)
 if(NOT BQPD)
    message(WARNING "Optional library BQPD was not found.")
 else()
-   list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp uno/solvers/BQPD/wdotd.f)
+   list(APPEND UNO_SOURCE_FILES uno/solvers/BQPD/BQPDSolver.cpp) # uno/solvers/BQPD/wdotd.f)
    list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/BQPDSolverTests.cpp)
    link_to_uno(bqpd ${BQPD})
 endif()

uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp

@@ -45,10 +45,6 @@ namespace uno {
 
    void QPSubproblem::evaluate_functions(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,
          const Multipliers& current_multipliers, const WarmstartInformation& warmstart_information) {
-      // Lagrangian Hessian
-      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
-         this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
-      }
       // objective gradient, constraints and constraint Jacobian
       if (warmstart_information.objective_changed) {
          problem.evaluate_objective_gradient(current_iterate, this->objective_gradient);
@@ -57,6 +53,10 @@ namespace uno {
          problem.evaluate_constraints(current_iterate, this->constraints);
          problem.evaluate_constraint_jacobian(current_iterate, this->constraint_jacobian);
       }
+      // Lagrangian Hessian
+      if (warmstart_information.objective_changed || warmstart_information.constraints_changed) {
+         this->hessian_model->evaluate(statistics, problem, current_iterate.primals, current_multipliers.constraints);
+      }
    }
 
    void QPSubproblem::solve(Statistics& statistics, const OptimizationProblem& problem, Iterate& current_iterate,  const Multipliers& current_multipliers,

uno/solvers/BQPD/BQPDSolver.cpp

@@ -14,14 +14,14 @@
 #include "options/Options.hpp"
 #include "fortran_interface.h"
 
-#define WSC FC_GLOBAL(wsc,WSC)
-#define KKTALPHAC FC_GLOBAL(kktalphac,KKTALPHAC)
-#define BQPD FC_GLOBAL(bqpd,BQPD)
+#define WSC FC_GLOBAL(wsc, WSC)
+#define ALPHAC FC_GLOBAL(alphac, ALPHAC)
+#define BQPD FC_GLOBAL(bqpd, BQPD)
+#define hessian_vector_product FC_GLOBAL(gdotx, GDOTX)
 
-namespace uno {
-#define BIG 1e30
+extern "C" {
+   void hessian_vector_product([[maybe_unused]] int *n, const double x[], const double ws[], const int lws[], double v[]);
 
-   extern "C" {
    // fortran common block used in bqpd/bqpd.f
    extern struct {
       int kk, ll, kkk, lll, mxws, mxlws;
@@ -30,13 +30,16 @@ namespace uno {
    // fortran common for inertia correction in wdotd
    extern struct {
       double alpha;
-   } KKTALPHAC;
+   } ALPHAC;
 
    extern void
    BQPD(const int* n, const int* m, int* k, int* kmax, double* a, int* la, double* x, double* bl, double* bu, double* f, double* fmin, double* g,
          double* r, double* w, double* e, int* ls, double* alp, int* lp, int* mlp, int* peq, double* ws, int* lws, const int* mode, int* ifail,
          int* info, int* iprint, int* nout);
-   }
+}
+
+namespace uno {
+   #define BIG 1e30
 
    // preallocate a bunch of stuff
    BQPDSolver::BQPDSolver(size_t number_variables, size_t number_constraints, size_t number_objective_gradient_nonzeros, size_t number_jacobian_nonzeros,
@@ -104,7 +107,7 @@ namespace uno {
       WSC.ll = static_cast<int>(this->size_hessian_sparsity);
       WSC.mxws = static_cast<int>(this->size_hessian_workspace);
       WSC.mxlws = static_cast<int>(this->size_hessian_sparsity_workspace);
-      KKTALPHAC.alpha = 0; // inertia control
+      ALPHAC.alpha = 0.; // inertia control
 
       if (this->print_subproblem) {
          DEBUG << "objective gradient: " << linear_objective;
@@ -311,4 +314,22 @@ namespace uno {
       }
       throw std::invalid_argument("The BQPD ifail is not consistent with the Uno status values");
    }
-} // namespace
+} // namespace
+
+void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]) {
+   for (size_t i = 0; i < *n; i++) {
+      v[i] = 0.;
+   }
+
+   int footer_start = lws[0];
+   for (int i = 0; i < *n; i++) {
+      for (int k = lws[footer_start + i]; k < lws[footer_start + i + 1]; k++) {
+         int j = lws[k] - 1;
+         v[i] += ws[k - 1] * x[j];
+         if (j != i) {
+            // off-diagonal term
+            v[j] += ws[k - 1] * x[i];
+         }
+      }
+   }
+}