From 1eed9ec8a18cf9fa022ac7d99c47d02e4a43a0fc Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 11 Dec 2024 19:17:34 +0100 Subject: [PATCH 1/2] BQPDSolver: moved the gdotx function from Fortran to C++ --- CMakeLists.txt | 2 +- uno/solvers/BQPD/BQPDSolver.cpp | 43 +++++++++++++++++++++++++-------- 2 files changed, 34 insertions(+), 11 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 716c45c2..7f5b3677 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) list(APPEND TESTS_UNO_SOURCE_FILES unotest/functional_tests/BQPDSolverTests.cpp) link_to_uno(bqpd ${BQPD}) endif() diff --git a/uno/solvers/BQPD/BQPDSolver.cpp b/uno/solvers/BQPD/BQPDSolver.cpp index 678780ea..c3ab4deb 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(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,24 @@ 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[]) { + assert(n != nullptr && "BQPDSolver::hessian_vector_product: the dimension n passed by pointer is NULL"); + + for (size_t i = 0; i < static_cast(*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 From 41f5d07b4e6712d9cc7746a9ac7e3a2ec19c1a4b Mon Sep 17 00:00:00 2001 From: Charlie Vanaret Date: Wed, 11 Dec 2024 19:20:03 +0100 Subject: [PATCH 2/2] Removed uno/solvers/BQPD/wdotd.f --- uno/solvers/BQPD/wdotd.f | 118 --------------------------------------- 1 file changed, 118 deletions(-) delete mode 100644 uno/solvers/BQPD/wdotd.f 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