Skip to content

Commit

Permalink
Merge pull request #135 from cvanaret/bqpd_gdotx
Browse files Browse the repository at this point in the history
BQPDSolver: moved the gdotx function from Fortran to C++
  • Loading branch information
cvanaret authored Dec 11, 2024
2 parents ec76b3a + 41f5d07 commit a68a05e
Show file tree
Hide file tree
Showing 3 changed files with 34 additions and 129 deletions.
2 changes: 1 addition & 1 deletion CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -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()
Expand Down
43 changes: 33 additions & 10 deletions uno/solvers/BQPD/BQPDSolver.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -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,
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -311,4 +314,24 @@ 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[]) {
assert(n != nullptr && "BQPDSolver::hessian_vector_product: the dimension n passed by pointer is NULL");

for (size_t i = 0; i < static_cast<size_t>(*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];
}
}
}
}
118 changes: 0 additions & 118 deletions uno/solvers/BQPD/wdotd.f

This file was deleted.

0 comments on commit a68a05e

Please sign in to comment.