diff --git a/CMakeLists.txt b/CMakeLists.txt index 716c45c2..f039b000 100644 --- a/CMakeLists.txt +++ b/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() diff --git a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp b/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp index fdba8a54..4c7e2511 100644 --- a/uno/ingredients/subproblems/inequality_constrained_methods/QPSubproblem.cpp +++ b/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, diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index 58aeea7a..ee3aceac 100644 --- a/uno/solvers/BQPD/BQPDSolver.cpp +++ b/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(this->size_hessian_sparsity); WSC.mxws = static_cast(this->size_hessian_workspace); WSC.mxlws = static_cast(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 \ No newline at end of file +} // 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]; + } + } + } +} \ No newline at end of file diff --git a/uno/solvers/BQPD/wdotd.f b/uno/solvers/BQPD/wdotd.f deleted file mode 100644 index e7ab5954..00000000 --- a/uno/solvers/BQPD/wdotd.f +++ /dev/null @@ -1,118 +0,0 @@ -C Copyright (c) 2018-2024 Sven Leyffer -C Licensed under the MIT license. See LICENSE file in the project directory for details. - -C hristen this file wdotd.f - - subroutine wdotd (n, x, ws, lws, result) - -c ========================================================== -c Computes result = W.x where W is Hessian and x is a vector for AMPL -c Assumes v=0 on entry (OK, if called from gdotx, see QPsolve*.f) -c ========================================================== - - implicit none - -c ... declaration of passed parameters - integer n, lws(0:*) - double precision x(n), result(n), ws(*) - -c ... declaration of internal variables - integer i, j, k, footer_start - -c inertia control for diagonal terms - double precision alpha - common /kktalphac/ alpha - -c ======================== procedure body ========================= - -c ... form result = W.x from sparse, upper triangular Hessian - footer_start = lws(0) - do i=1,n - do k=lws(footer_start + i - 1), lws(footer_start + i)-1 - j = lws(k) - result(i) = result(i) + ws(k)*x(j) - if (j.ne.i) then -c off-diagonal term - result(j) = result(j) + ws(k)*x(i) - else -c diagonal term - result(i) = result(i) + alpha*x(j) - endif - enddo - enddo - return - end - -c ****************************************************************** - - subroutine gdotx (n, x, ws, lws, result) - - implicit none - -c ... declaration of passed parameters - integer n, lws(*) - double precision x(n), result(n), ws(*) - -c ... declaration of internal variables - integer i - -c ... storage map for hessian and scale_mode - integer scale_mode, phe - common /scalec/ scale_mode, phe - -c ======================== procedure body ========================= - -c ... set result = 0 - do i=1,n - result(i) = 0.D0 - enddo - -c ... allow for scaling of variables - if ((scale_mode.eq.1).or.(scale_mode.eq.3)) then - do i=1,n - x(i) = x(i) * ws(i) - enddo - endif - -c ... form v = W.d from sparse, upper triangular Hessian - call Wdotd (n, x, ws(phe+1), lws, result) - -c ... allow for scaling of variables - if ((scale_mode.eq.1).or.(scale_mode.eq.3)) then - do i=1,n - result(i) = result(i) * ws(i) - x(i) = x(i) / ws(i) - enddo - endif - - return - end -c ****************************************************************** - subroutine saipy2(s,a,la,i,y,n) - implicit double precision (a-h,o-z) - dimension a(*),la(0:*),y(*) -c ======================== procedure body ========================= -c saxpy with column i of A: y + s*A_{i, :} - if(s.eq.0.D0) return - j_column_start = la(0) + i - do j = la(j_column_start), la(j_column_start+1)-1 - i_variable = la(j) - y(i_variable) = y(i_variable) + s*a(j) - enddo - return - end - -c **************************** E N D ********************************* - function daiscpr2(n,a,la,i,x,b) - implicit double precision (a-h,o-z) - dimension a(*),la(0:*),x(*) - DOUBLE PRECISION daiscpr2 -c dot product of x and row i of A - daiscpr2 = dble(b) - j_column_start = la(0) + i - do j = la(j_column_start), la(j_column_start+1)-1 - i_variable = la(j) - daiscpr2 = daiscpr2 + dble(x(i_variable))*dble(a(j)) - enddo - return - end