Moved gdotx function (Hessian-vector product) from Fortran file to BQ…

…PDSolver.cpp

CMakeLists.txt

@@ -129,7 +129,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()

uno/solvers/BQPD/BQPDSolver.cpp

@@ -17,8 +17,11 @@
 #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)
 
 extern "C" {
+   void hessian_vector_product(int *n, const double x[], const double ws[], const int lws[], double v[]);
+
    // fortran common block used in bqpd/bqpd.f
    extern struct {
       int kk, ll, kkk, lll, mxws, mxlws;
@@ -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 (int 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];
+         }
+      }
+   }
+}